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COMPUTER PROGRAM FOR STATIC AND DYNAMIC AXISYMMETRIC 
NONLINEAR RESPONSE OF SYMMETRICALLY LOADED 
ORTHOTROPIC SHELLS OF REVOLUTION 

By Wendell B. Stephens 
Langley Research Center 

SUMMARY 

A computer program has been developed which determines the nonlinear behavior 
of symmetrically loaded elastic orthotropic shells of revolution. The loading can be due 
to either mechanical or thermal forces and can be applied statically or dynamically. The 
analysis is based on Sanders' equations for shells with small strains and moderately 
small rotations and allows for variable stiffness properties of the shell along the merid- 
ian. Spatial derivatives are approximated by finite differences, and integration with 
respect to time is carried out by the Houbolt method. For static behavior, or dynamic 
response at each point in time, a Newton-Raphson method is applied for convergence to 
the nonlinear solution. The boundary conditions are presented in a general form which 
allows either classical or elastic constraints to be used. The program, which is written 
in FORTRAN IV language, is described in detail and sample calculations are included. 

INTRODUCTION 

The analysis of shells of revolution subjected to static, thermal, or time-dependent 
loads is an important problem in the design of missiles and space vehicles. A finite- 
difference solution for the linear bending behavior of an isotropic shell subjected to an 
arbitrary static load is contained in reference 1 and is modeled after the analysis pro- 
cedure found in reference 2. Geometrically nonlinear terms are included for essentially 
the same problem in reference 3. However, there remains a need for a program which 
accounts for dynamic loads and material orthotropy. Such a dynamic response analysis 
is useful for practical aerospace applications such as the study of launch, staging, and 
water-impact loadings of aeroshells. In addition, such an analysis would provide a means 
of determining the nonlinear prebuckling stress distributions required for accurate sta- 
bility analyses. In this report a computer program is described which has been developed 
to determine the axisymmetric nonlinear static and dynamic response including axisym- 
metric static and dynamic buckling of an arbitrary elastic orthotropic shell of revolution 



subjected to axisymmetric loads. The analysis, programing techniques, and the computer 
program documentation are presented as well as representative sample problems. 

The analysis is based on Sanders' nonlinear equations (ref. 4) with material orthot- 
ropy added as in reference 5. The governing partial differential equations are written in 
terms of first-order spatial derivatives and solved numerically by using central differ- 
ences for derivatives along the meridian and backward differences for time derivatives. 
Integration with respect to time is started by using the Houbolt technique (refs. 6 and 7). 
For the boundary conditions, either classical or elastic constraints may be used. The 
nonlinear difference equations are solved for each time step or static load increment by 
the Newton-Raphson method (ref. 8). "Top-of-the-knee" static buckling is determined 
from the lack of convergence of the Newton-Raphson procedure. 

The program is divided into nine subroutines and seven user-supplied function sub- 
programs. A maximum of 101 equal stations is provided requiring an octal storage of 
70 000 memory words. The program is written in CDC version of FORTRAN IV language 
for operation in the CDC 6600/6400 digital computer at the Langley Research Center. 

The output consists of a problem description together with displacements, rotations, and 
moment and force resultants in tabular form. 

In order to present both the analysis and the computer program, appendixes are 
frequently used to simplify the text. Appendixes A, B, C, D, and E are used to clarify 
the presentation of the analysis, and appendixes F and G contain the program listing and 
a sample of the program output, respectively. 

SYMBOLS 


The units for physical quantities defined in this paper are given both in the U.S. 
Customary Units and in the International System of Units (SI). Factors relating the two 
systems are given in appendix H. 

a reference (or characteristic) length 


C 11 ’ C 12> C 22 


nondimensional orthotropic extensional material constants defined in 

Cn 


equations (14); for example, C-q = 


EqHq 


C 11 ’ C 12> C 22 

D 11 ’ D 12’ D 22 


orthotropic extensional material constants 


nondimensional orthotropic bending material constants defined in 

>! ®ii 


equations (14); for example, D-q = — ^ 


2 



D 11 ’ D 12> D 22 


orthotropic bending material constants 


E 0 reference modulus of elasticity 

Ej,E 2 moduli of elasticity in principal directions, meridional and circumferential, 
respectively 

ElO’^20 nondimensional moduli of elasticity in principal directions; for example, 

Ei 

*10 “l 1 

“O 

E 0 /pg specific stiffness where pg is weight density 

e ll’ e 22 nondimensional extensional strain, meridional and circumferential directions, 
respectively; for example, e-Q = e ii>l 

g acceleration due to gravity 

H shell thickness 


H 0 reference thickness 

H maximum shell rise of spherical cap considered in sample problem 

h nondimensional shell thickness, H/Hq 


j 


temperature exponent in equation (17) 


K 11 ’ K 12’ K 22 


nondimensional orthotropic material constants associated with 
coupling between extension and bending and defined in equa- 


tions (14); for example, K-jj 


xK n 


Kh,Ki 2,K22 orthotropic material constants associated with coupling between 

extension and bending 


^ 1 1*^22 principal change in curvatures, meridional and circumferential directions, 
respectively 
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Mh,M 22 bending-moment resultants in principal directions, meridional and circum- 
ferential, respectively 

m ll> m 22 nondimensional bending-moment resultants in principal directions; for 

aMii 

example, m^j = — 

oiii 


Nn,N22 membrane stress resultants, meridional and circumferential directions, 
respectively 


n 


number of stations along meridian 


n H , n 22 nondimensional membrane stress resultants, meridional and circumferential, 


respectively; for example, njj = 


N 
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crHo 


P,P S 


lateral and meridional forces per unit area, respectively 


cr 


critical symmetric buckling load 


p* = _E_ 
p cl 


p,p s nondimensional lateral and meridional forces per unit area, respectively; 

Pa 


for example, p = 


ctH r 


p cl 


nondimensional classical buckling pressure of complete spherical shell 
(see eq. (30)) 


Q transverse shear resultant 

q nondimensional transverse shear resultant, Q/ aHo 

R radial distance from axis of symmetry to shell reference surface 

Rj,R 2 principal radii of curvature, meridional and circumferential directions, 

respectively 
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r 


r l’ r 2 


nondimensional radial distance from axis of symmetry to shell reference 
surface, R/a 


nondimensional principal radii of curvature; for example, r j 


Rl 


a 


S 

s 

As 

t,t 1 ,t 2 

t 

At 

,m ,m 

tl ,t 2 


U,W 

u,w 


X 

y 

z 




a,/3,y,6 


distance measured along shell meridian 
nondimensional distance measured along meridian, S/a 
nondimensional meridional difference increment 
temperature quantities defined with equation (17) 
real time 

real time increment 

nondimensional thermal moment resultant in principal directions, defined in 
equations (16) 

nondimensional thermal force resultant in principal directions, defined 
in equations (16) 

meridional and normal displacement, respectively 

nondimensional meridional and normal displacement, respectively; for 
example, u = ^ 

force vector with elements njj, q, and m^ 
displacement vector with elements u, w, and /3 
vector composed of x and y vectors 

coefficients of linear thermal expansion in principal directions, meridional 
and circumferential, respectively 

coefficients of the acceleration difference equation (21) 
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0 

0 

A 

6 pg 

e ll ,e 22 

? 

V 

e 

* 11’ *22 

\ 

^ 12^21 

P 

a 

T 

At 

Subscripts: 


nondimensional rotation, ?7/3 


meridional rotation 

average deflection of spherical cap in sample problem 
Krone cker delta 
membrane strains 

coordinate normal to reference surface of shell, positive outward, with origin 
at reference surface and nondimensionalized by H 0 

ratio of reference elasticity modulus to reference stress, E 0 /a 

circumferential coordinate 

nondimensional principal curvatures; for example, Kjj = a^k^ 
ratio of reference thickness to characteristic length, H 0 /a 
shell parameter defined by equation (31) 

Poisson’s ratios for meridional and circumferential directions, respectively 
mass density 
reference stress 


nondimensional time, t 



nondimensional time increment 


colatitude angle, angle between shell axis and normal to shell middle surface 


spatial station number, that is, 1, 2, . . ., n 
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j matrix number, that is, 1, 2, . . 2n 

k kth equation of set of equations at a point 

m time step, that is, 1, 2, . . . 

max maximum 


Matrices: 


z,e 6x1 


A,B,C,D,E7i 
F1jF 3 ,P,Q J 


3X3 


XOLD,X,"l 
R,x,y,q,Z J 


3X1 


H,H,M 6x6 


A prime indicates a derivative with respect to the nondimensional meridional 
distance s. 

A dot indicates a derivative with respect to nondimensional time t. 


ANALYTICAL FORMULATION 


The shell analysis procedure is summarized in this section. Also included are the 
geometric description, derivation of the nonlinear equilibrium conditions, compatibility 
equations, and differencing scheme as well as the Newton-Raphson procedure for solution 
of the governing equations. 


Shell Geometry 

The shell geometry and coordinate system for the reference surface of a general 
shell of revolution are shown in figure 1. The geometry of the shell reference surface is 
defined by <b and R. Any point in the shell may be located by specifying the orthogonal 

a 

coordinates s, 8, and £ where s = ^ and is the nondimensional meridional coordinate, 

a 

S is the meridional shell coordinate, a is the reference length of the shell, 6 is the 
circumferential coordinate, and £ is a coordinate normal to and originating at the shell 
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reference surface, positive outward. The nondimensional principal radii of curvature, 
rj and r 2 , are (ref. 1) 



1 _ sin 0 

r 2 " r J 


( 1 ) 


where the prime indicates a derivative with respect to the nondimensional meridional 
distance s. The radii are nondimensionalized by use of the reference length a. 


Equilibrium Conditions 

By utilizing the results in reference 4 and the nondimensional variables described 
in reference 2, the nondimensional equilibrium equations become 


cos <2> cos 0 1 , fn 

n ll + n H + ^ q p— n 22 + - <t> 'fan - hu = -P s 


cos cp sin <p 


q' - <p'n n + — ~ q - 


r n 22+^ 


fall) + y- /3n n J - hw = -p 


COS (j) 

'it + ~r 


COS <p 


m,, + — — m n — m 22 - ~s = 0 


( 2 ) 

(3) 

(4) 


where njj and n 22 are the stress resultants, mu and m22 are the bending- 
moment resultants, u and w are the displacements, p s and p are the surface 
loads, q is the shear resultant, and /3 is the meridional rotation. These quantities 
are defined in figure 2. The term X is a nondimensional constant representing the 
ratio of the reference thickness Ho to the reference length a. The dots indicate deriv- 
atives with respect to nondimensional time t. 

Rotational and Strain-Displacement Relationships 

The nondimensional rotation, strain-displacement relationships, and curvatures 
from reference 4 are 


/3 = w' - (p’u 


(5) 


e n = u’ + 0-W+J-/3 2 


( 6 ) 
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cos 4> sin <P 

e 22 = u + r~ w 


( 7 ) 


“11 = -i 3 ' 

cos <b 
*22 = F~ 


/3 


(8) 

(9) 


The terms e-^ and e 22 are the nondimensional principal strains and and K 22 

are the nondimensional principal curvatures. 

Constitutive Equations 

For symmetrically loaded orthotropic shells of revolution nondimensional elasticity 
relationships obtained from reference 5 can be written as 


n 


n ll _ C ll e ll + c 12 e 22 + K ll /C n + K 12*22 " *1 


( 10 ) 


n 


n 22 “ C 12 e ll + C 22 e 22 + K 12*ll + K 22 k 22 “ *2 


(ID 


m H - ''•”^l^ll e ll + ^"^1^12 e 22 + *~^ll K ll + *”^12*22 " ^ (12) 


m 22 ~ X ~ 2K 12 e ll + X ~ 2K 22 e 22 + X ~ 2d 12 K H + x " 2d 22 k 22 " 1™ ( 13 ) 


Since only axisymmetric behavior is considered, only these four relationships are 
required. The nondimensional stiffnesses are given by 

■\ 


C 


11 


E l° - f ?2 d C 
1 - v 12 v 21 "?1 


c 12 


c 22 


Kll 


^i 2 E io r? 2 d? 

1 - V21 j ?i 

Eg ° d C 

1 - ^12^21 5l 

AE 10 p?2 

— \ 5 d? 

1 ” ^12^21 ? 1 


(14) 


(Equations continued on next page) 
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Xi/ 12 E 10 

k 12 „ 

r? 2 

C d? 

1 “ v V2 v 2 1 

J ?i 

„ aE 20 

r? 2 

k 22 

? d? 

1 " v 12 v 21 

J ?1 

x2e 10 

D 11 - „ 

f Ca ? 2 d? 

1 " ^12^21 

J ?1 

X 2 i^i2 E io 

d 12 = 

f 2 ? 2 d? 

1 ‘ ^12^21 ' 

J ?1 

A 2 E 2 q 

D 22 = — 

r ?2 £>2 di * 

1 - V 12 1J 21 ' 

j ?i 




( 14 ) 


where £ is positive outward and and ? 2 are the distances to the inner and outer 
shell surfaces, respectively, from the reference surface. 

Because of the symmetry of the orthotropic constants, use has been made of the 
relationship 


e 10 p 12 = e 20^21 


( 15 ) 


The nondimensional thermal forces and moments in the meridional and circumferential 
directions, respectively, due to a temperature T(s,£) are (ref. 1) 



E 10 ? l 

1 - ^12^21 



T d? 




,n 

E 20 t ? / 

\ r ? 2 

T d£ 

i 2 ~~r 

(« 2 

+ "2i“i L 

i 

- "12 "21' 

' J Ci 


t m -- 
1 i 

El °” h 

" "l2"21 V 

+ 

T? d£ 

t m - 

2 T 

E 20 7 l /„, 

( a 2 

~ u 12 "21' 


n d? 


( 16 ) 
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where a j and are the orthrotropic coefficients of thermal expansion in the prin- 
cipal directions. 


Temperature Profile 

The temperature is allowed to vary through the thickness of the shell and along the 
meridian as follows 

T(s) = T 1 (s) + T 2 (sK j (17) 

where the Tj defines the temperature change from a standard temperature at the ref- 
erence surface and T2 is the difference between the temperatures of the shell outer 
and inner surfaces at ? 2 and ?l- The exponent j is used to define the temperature 
thickness profile as a constant (j = 0) or as a linear variation through the thickness (j = 1) 
or as a nonlinear variation through the thickness (j =2). 

Finite -Difference Formulation of Governing Equations 

It is shown in appendix A that equations (2) to (13) can be written as six partial dif- 
ferential equations. These equations are first order in spatial derivatives and second 
order in time derivatives. The set of equations in matrix form is 

Iz' + (h + H)z = e + Mz (18) 

where 


r 

n ll 

q 

, m n v 
z = < 

u 

w 

13 

v J 

Here z is the solution vector of six variables, I is the 6x6 identity matrix, H and 
H are the linear and nonlinear 6x6 coefficient matrices of z, respectively, M is the 
6x6 mass matrix of z, and e is the six-element load vector. The elements of H, H, 
e, and M, are listed in appendix A. 

The governing equations are converted into difference equations by utilizing central 
differences for the spatial derivatives and backward differences (refs. 6 and 7) for the 
time derivatives. As shown in reference 7, this backward-difference scheme is 
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numerically stable. The spatial finite-difference representations are written at a point 
halfway between stations as shown in figure 3 and are of the form 


z i-l/2 


z i + z i-l 


(19) 


z i - 1/2 


z i - z i-l 
As 


( 20 ) 


The second-order time derivative in equations (18) is approximated at the ith station by 


z i,m 


(At) 


2 ( amZ i,m + ^m z i,m-l + y m z i,m-2 + 6 m z i )m -3 


( 21 ) 


where i = 1, 2, . . n and m = 1, 2, . . . . In equations (19), (20), and (21) the sub- 
scripts i and m indicate spatial and time stations, respectively, and As and At 
are the spatial and time increments, respectively. The coefficients 'a m , 0 m , y m , and 
6 m depend on the time step and initial conditions and are given in appendix B. Applica- 
tion of these finite-difference approximations (eqs. (19) to (21)) to the governing equa- 
tions (18) leads to the following set of nonlinear algebraic equations at the mth time step: 


F i-l/2 z i-l,m + G i-l/2 z i,m ~ L i-l/2 


where 


F i-l/2 
G i - 1/2 

L i - 1/ 2 


1 

2 




+ ^i - 1/2 

+ fi i-l/2 


M i-l/2%\ I 
(At ) 2 ) As 

M i-l/2%\ + J_ 
(AT ) 2 ) As 






^i-1/2/— — 

e i- 1/2 + (a ’ t) 2 v mZ i-l/2,m-l + y m z i-l/2,m-2 

+ ^m z i-i/2, m -3) 


(22) 


(23) 


and i = 2, 3, . . ., n and m = 1, 2, 3, . . . . These equations with three boundary con- 
ditions at each edge of the shell define the problem to be solved and must be solved 
simultaneously to determine z at the mth time step. 
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Boundary Conditions 


As shown in reference 4 the classical shell boundary conditions at either edge, 

g 

s = 0 and s = ma ~^ , are defined by either force resultants njj, q, and mu or dis- 

cL 

placements u, w, and /3, so that 


f2xj + Ay i = l 


(24) 


Here the subscript i is either 1 or n and the 3x3 O and A matrices and the 
3x1 l vector are required to define the boundary conditions. The vectors Xj and 
yi are 


r n ir 


f U ] 

q 

<< 

11 

w 


i 

i 



(25) 


These vectors define the force and displacement subvectors of z, respectively. Typical 
boundary conditions including general elastic constraints are discussed in appendix C. 


Computational Procedure 

The nonlinear set of equations (22) and (24) are linearized by use of an iterative 
Newton-Raphson procedure (ref. 8). This is done by placing the Li_i/2 term and the 
l term on the left-hand side of equations (22) and (24), respectively, and writing the kth 
equation at the ith station as 

f k( z i> z i-l> s ) = 0 (26) 

where k is 1, 2, . . ., 6 for equation (22) and k = 1, 2, 3 for equations (24). Use of the 
first two terms of a Taylor's expansion for equation (26) together with an approximate 
solution vector z gives 


f k( z i> z i_i> s ) = f k( z i> z i_i> s ) 


d h 

6Zi + — ! 

1 a z , 


z i= z i 


i-1 


6z i-l = 0 


4-l =z i-l 


(27) 


where i = 1, 2, . . ., n and where 6z^ is the correction vector which must be added to 
the approximate solution vector zj so that equation (26) is satisfied. 
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The iterative procedure consists of adding the correction vector to the approximate 
solution vector to obtain an improved approximate solution. Thus 

zj +1 = zj + 6zj (28) 

where the superscript j indicates the jth iteration cycle. When 5z^ becomes suffi- 
ciently small, convergence has been obtained. 

For convenience in the solution of the simultaneous equations, the correction vector 
6zj is partitioned into the two three-element ordered subvectors 5x^ and 6yj. Thus 
the set of governing equations (27) including the appropriate boundary conditions take the 
following form of a five-diagonal-banded matrix where each element is a 3 x 3 matrix. 


r* 

Cl 

Dl 

Ei 










^1 

b 2 

c 2 

d 2 

e 2 







yi 


V-2 

As 

b 3 

c 3 

d 3 

e 3 






x 2 


q 3 



A 2(i-1) 

B 2 (i-1) 

c 2 (i— 1) 

D 2(i-1) 

e 2 (i— 1) 



< 

x i 

> = < 

l2(i-l)> 




A 2i-1 

B 2 i-1 

C 2i -l 

D 2i-1 

E 2i-1 



yi 


021-1 






A 2n-1 

B 2n -1 

C 2n-1 

D 2n-1 


x n 


02n-l 







A 2n 

B 2n 

c 2n 


y n 
^ J 




(29) 

For brevity, 6xj and 6yj are written as x^ and y^ in equation (29). The first 

g 

and last rows are the boundary conditions at s = 0 and s = respectively, and are 

obtained from equations (24). The six first-order governing equations from equations (18) 
and (27) correspond to the pair of rows at 2(i-l) and 2i-l, respectively. Here n is 
equal to the number of spatial stations. The definitions of the A, B, C, D, E, and 
q matrices in terms of equations (2) to (13), (22), and (27) are given in appendix D. 
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The set of equations (29) is solved by using a modified Potters method (refs. 2 or 9) 
for banded matrices. A presentation of the recurrence equations required for the Potters 
method is contained in appendix E. For each time step, the elements in the 3X3 matri- 
ces A, B, C, D, and E and the three-element vector q are functions of the shell 
properties, and the new displacement and stress state for the last three time steps. 

The vectors z^q and Zio at the initial time t = 0 must be given. Both the 
incorporation of the initial conditions into the problem and the definition of the 'a m , /3 m , 
y , 6 m coefficients for z m are contained in reference 7 and appendix B. Appen- 
dix B also includes nonhomogeneous initial conditions. 

This analysis and numerical solution has been programed in FORTRAN IV and the 
resulting computer program (SADA0S) includes the input provisions of general shell 
shape, thermal and mechanical loads, structural orthotropy, and arbitrary boundary con- 
ditions at each end of the shell. 


COMPUTER PROGRAM 

This section contains the description of the computer program SADA0S and is 
intended to be a user's document. A listing of the program is contained in appendix F 
and a sample printout in appendix G. In writing this program, various options on types 
of analyses, geometry, and boundary conditions have been included to eliminate the nec- 
essity of having the user develop subroutines. However, should these options be inade- 
quate, the program is subdivided into separate subroutines so that further options can be 
exercised without a detailed knowledge of the program. Certain function subprograms 
are required to be programed by the user. These function subprograms define the 
loading, shell thickness, temperature integrals, and initial conditions. In addition, input 
data and computer subroutine preparation are explained in detail in later sections. 

Program Organization 

The flow chart is presented in the following block diagram. As an aid in reading 
the block diagram, a list of subroutines and their description is presented in table 1. In 
table 2 the variables and constants are listed with their program names. 
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( 1 ) 

( 2 ) 


Read Namelists 

I 


Call GE0MTY 


(3) 




Call B0C0N 1 


(4) I Call ABCDES 


If IDYM = 1 
Call DYNAMIC 


if IDYM = 2 and no 
convergence after 20 
attempts then P = AP/5| 


(5) 

( 6 ) 
(7) 


Call STIF 1 


Call FUNCT 


P = P + AP 


[Call POTTER 
{Convergence) 



( tnd of static 
analysis 


(9a) (8) 1 Call ANSWERS 


(9c) 


^V ftDYM = O) — (IDYM> — ( IDYM = 2 

{9b) 4 dym = l) 


Has P been decrease 
times 


Calculate z and z 




K0UNT = K0UNT + 1 


- <Is K0UNT = KMAX) 



End of static' 
buckling > 


End of dynamic 
.analysis 


Block diagram of SADA0S 
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A detailed description of the computing in each block of the block diagram is as follows: 


BLOCK (1): 

BLOCK (2): 
BLOCK (3): 
BLOCK (4): 


BLOCK (5): 


BLOCK (6): 


BLOCK (7): 


BLOCK (8): 
BLOCK (9a): 
(9b): 


Namelist GIVEN containing basic input data is read in. If requested, the 
optional information on boundary conditions will be read in through name- 
lists ELB0L and ELB0R. 

Shell geometry is defined at the i-1/2 increment midpoints along the shell 
meridian. The values defined by GE0MTY are r, <p, and <p'. 


Boundary conditions at each end are set. Matrices Cj, Dj, Ej, q-^, 
A 2n> B 2n> c 2n> and <l2n in equation (29) are defined. 

Matrices A, B, C, D, E, and q from equation (29) are calculated. 
These matrices are further defined in appendix D. These matrices are 
calculated for each i = 2, 3, . . ., n and, in turn, call the subroutines 
DYNAMIC (if IDYM = 1) and STIF. STIF sets the C, D, and K values 
from equations (14). DYNAMIC sets a, j§, y", 6 of equations (21). 

From equation (27) the term fk(%>^i-l> s ) is calculated and placed on the 
right-hand side of that equation. This corresponds to the q vector in 
equation (29). 


The Potters method or Gaussian elimination scheme described in appendix E 


is used to solve equation (29) for the vector 6z where 6z 
improved approximate vector is given by equation (28). 


f ® x i\ 
1 ~ WiJ 


The 


If the norm of 6z is very small compared with the norm of the improved 
vector z, then the problem has converged. If 6z is not sufficiently small, 
a return is made to block (3), and blocks (3) to (7) are repeated until con- 
vergence is obtained. If a static buckling problem (IDYM = 2) is being 
attempted, and after 20 iterations there is still no convergence, a return is 
made to the last converged load solution and a smaller load increment step 
is attempted. 


After convergence there is a tabular printout of the following variables at 
each station i: n^, q, mji, u, w, /3, n22, m 2 2 > P, and p s . 


If the analysis is a static stress analysis (IDYM = 0), the program 
terminates. 


If the analysis is a dynamic response (IDYM = 1) problem, the z and z 
vectors are calculated and an increment in time and time step is taken. 
This procedure is continued until KMAX steps are taken. 
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(9c): If static buckling load (IDYM = 2) is desired, the load is increased. If 

from block (7) the load increment has not been decreased five times, the 
problem returns to block (3). If five load increment reductions have taken 
place, the problem terminates. 


Input Data 

The following quantities in namelist GIVEN must be defined: N, NTYPE, CHAR, 
H0, S0, R0, PHI0, PP, PPS, E0, El, E2, NU12, NU21, SIG0, N0NLIN, C0NV, IDYM, 

KM AX, DTAU, ALFA1, ALFA2, Tl, T2, ITEMP, LBCL, LBCR, SL, SR, IFREQ, and 
ISTART. The format for the input data contained in a namelist is given in reference 10. 
The first column of the data cards cannot be used. The definitions of these quantities 
are as follows: 

Name Type Interpretation 

N integer n, number of stations along the meridian 

NTYPE integer sets type of shell geometry to be analyzed: 

NTYPE = 1 denotes a cylindrical shell. 

R0 is the radius from the shell axis to the shell reference 
surface. 

PHI0, the colatitude angle, is 90°. 

S0 is the shell length. 

NTYPE = 2 denotes a conical shell. 

R0 is the radial distance to the first station at S = 0 from the 
shell axis. 

S0 is the length along the shell meridian. 

PHI0 is the colatitude angle (i.e., 90° — Semivertex angle). 
NTYPE = 3 denotes a spherical cap. 

R0 is the shell radius. 

PHI0 is the colatitude angle at S = S max . 

S0 is calculated internally and is read in as zero. 

NTYPE = 4 denotes that the user will read in a special geometry 
by adding statements to GE0MTY as required to define r, <p, and 
cj>' at each i - 1/2 station. Input constants R0, PHI0, and S0 can 
be used as desired by the programer. The statements are placed 
after the card labeled 50 and before the card labeled 60. 


CHAR 

real 

a, reference shell dimension to be selected by user and used inter 
nally for nondimensionalizing the geometry and output quantities 

H0 

real 

H 0 , reference thickness 
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Name 


S0 

R0 

PHI0 

PP 

PPS 

E0 

El 

E2 

NU12 

NU21 

SIG0 

N0NLIN 

C0NV 


IDYM 

KMAX 


DTAU 

ALFA1 

ALFA2 


Type 

real 

real 

real 

real 

real 

real 

real 

real 

real 

real 

real 

integer 

real 


integer 


integer 


real 


Interpretation 

input quantity defined by NTYPE 
input quantity defined by NTYPE 

input quantity defined by NTYPE and read in degrees 

constant used to define p, the normal pressure, in FUNCTION 
PL(I) 


constant used to define p s , the meridional pressure, in 
FUNCTION PS(I) 

E 0 , reference elasticity modulus 

Ej_, elasticity modulus in meridional direction 

E 2 , elasticity modulus in circumferential direction 

^ 12 ’ P°i sson ' s ratio in the meridional direction 

v 21> Poisson's ratio in the circumferential direction 

a, reference stress level; normally SIG0 = 1. 

If a linear solution is desired, set N0NLIN = 0. If nonlinear 
terms are to be included, set N0NLIN = 1. 

convergence criteria. Compares the error norm with the norm of 
the approximate solution vector (i.e., ), to insure conver- 

\ II z 11 / 

gence to proper order of magnitude. Usually C0NV = 1 . x 10“^. 

If static stress analysis is desired, set IDYM = 0. If dynamic 
response analysis is desired, set IDYM = 1. If static buckling 
analysis is desired, set IDYM = 2. 


number of time steps desired when IDYM = 1. Maximum number 
of static-load solutions when IDYM = 2. (Provides an upper limit 
on iterations when snap-through buckling (IDYM = 2) does not 
occur.) 


At, size of the nondimensional time increment. 



real oij, coefficient of thermal expansion in the meridional direction 

real 05 , coefficient of thermal expansion in the circumferential 

direction 
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Name 


Interpretation 


Type 

T1 real constant used to define Tj in equation (17) 

T2 real constant used to define T 2 in equation (,17) 

ITEMP integer j, integer exponent used in equation (17) 

LBCL integer sets boundary condition at the i = 1 (S = 0) edge (see eqs. (C2) 

to (C5)): 

LBCL = 1 is a pole point 
LBCL = 2 is a pinned edge 
LBCL = 3 is a fixed edge 
LBCL = 4 is a free edge 

LBCL = 5 elastic constants in namelist ELB0L must be given 

LBCR integer sets boundary conditions at the i = n (S = S max ) edge (see 

eqs. (C2) to (C 5)) : 

LBCR =1 is a pole point 
LBCR = 2 is a pinned edge 
LBCR = 3 is a fixed edge 
LBCR = 4 is a free edge 

LBCR = 5 elastic constants in namelist ELB0R must be given 

SL real three-element array equated to values defined by LBCL (see 

appendix C), equivalent to Zj in equations (24) and (Cl) 

SR real three-element array equated to values defined by LBCR (see 

appendix C), equivalent to Z n in equations (24) and (Cl) 

IFREQ integer frequency of printout at time steps of dynamic-load problems 

(IDYM = 1) or at load steps in the static buckling problem 
(IDYM = 2) 

ISTART integer Normally ISTART = 0. If ISTART = 1, user must supply non- 

homogeneous initial values to FUNCTION DV for deflections and 
velocities u, w, u, and w. 

The first input quantity in namelist GIVEN is preceded by $GIVEN and the last 
input quantity is followed by $. For example, the first input card could be 

$GIVEN N = 21, NTYPE = 3, H0 = 1., 

and the last namelist card could be 

SL(1) = 0., 0., 0., SR(1) = 0., 0., 0., IFREQ = 4, ISTART =0$ 
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Finally, one input card may contain a description of the problem. All 80 columns 
may be used. If no description is desired, a blank card must be included after the last 
namelist GIVEN data card. 

If LBCL = 5, then namelist ELB0L must be included in the input. Two 3x3 matri- 
ces Gl and Al are read in columnwise. Their elements specify the elastic con- 
straints at the first boundary of the meridian (i = 1). For example, a simply supported 
edge of a shell free to displace in a horizontal plane and with an applied edge moment 
yields the boundary conditions 

nu cos 0 O + 1 sin 0o = 0 
u sin <p o - w cos <p 0 = 0 

“11 = m 0 

where <fr Q is the colatitude angle at the boundary. Thus, the S2, A, and SL matrices 
become 



If LBCR = 5 the namelist ELB0R must be included in the input. Two 3x3 matri- 
ces fin and Ar are read in columnwise. Their elements specify the elastic con- 
straint at the last boundary of the meridian (i = n). 

User-Prepared Function Subprograms 

In addition to the input data, the user must prepare certain function subprograms. 
These function subprograms must be written and included in the program by the user to 
calculate the quantity at each half-station, i-1/2. For example, at i = 2 the FUNCTION 
PL(I) will define the lateral nondimensional pressure at a point halfway between i = 1 
and i = 2. The i-1/2 station is shown in figure 3. The following table describes the 
function subprograms. 

FUNCTION Quantity Interpretation 

PL(I) Pi-1/2 (i = 2, 3, . . ., n) computes the nondimensional lateral 

pressure 
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FUNCTION 

Quantity 

Interpretation 

PS(I) 

p s,i-l/2 = 2 > 3 > • • •» n ) 

computes the nondimensional meridional 
pressure 

TDZ (I) 

J ? ?2T i-l/2 d ? ft = 2, 3, . . ., n) 

computes the nondimensional thermal 
force integral 

TZDZ 

|a 

2 T i_l/ 2 ? d ? ft = 2, 3, . . ., n) 

computes the nondimensional thermal 
moment integral 

T(I) 

h i-l/2 ft = 2, 3, . . ., n) 

computes the nondimensional shell 
thickness 

IP(K0UNT) 

required when IDYM = 1 

locates time stations where a load is sud 
denly applied or removed 

DV(M,I) 

required when 1ST ART = 1 

prescribes the nonhomogeneous initial 
conditions: 


M = 1 denotes the u displacement at 
station i 

M = 2 denotes the w displacement at 
station i 

M = 3 denotes the u velocity at 
station i 

M = 4 denotes the w velocity at 
station i 

Program Output 

The output is divided into two parts. The first part is a printout of the input data 
and shell geometry. The second part is the printout of njj, q, mjj, u, w, /3, p, 
p S) n22, and m22 at all stations for the converged solution. If the problem varies with 
time, then the second part is repeated KMAX times. At a time station where there is a 
sudden change in load (i.e., step loads denoted by IP = 1 in FUNCTION IP(K0UNT)) 
there is an additional printout of the vectors u, w, u, and w. 

Program Limitations 

The program is limited to 101 spatial stations and 70 000 octal storage locations. 

At present there is no programed mechanism for allowing the orthotropic coefficients of 
thermal expansion aj and a*j to vary through the thickness. This could be accom- 
plished by the user by writing FUNCTION TDZ and FUNCTION TZDZ to include variable 
thermal coefficients ai(?) and a*>(£). 
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In subroutine STIF these constants are set up for a general orthotropic shell with 
the reference surface at the middle surface. The stiffnesses, that is, the Cij, Dij, and 
Kjj defined in equations (14), are general. Thus, the user could, with minimal knowledge 
of the computer program and equations (14), alter subroutine STIF to include a shell stiff- 
ened by rings or stringers smeared over appropriate increments (ref. 11) and using any 
reference surface. 


Analytical Limitations 

For static buckling problems (IDYM = 2) the buckling is limited to "top-of-the-knee" 
axisymmetric buckling loads. A detailed discussion of top-of-the-knee buckling is con- 
tained in reference 12. 

Errors in results will normally be one of two types: (1) inconsistency of input data 

or (2) numerical error inherent in the analysis. The first type of error can be eliminated 
by careful scrutiny and checking of the input data and user-prepared function statements. 
The second type of error can only be minimized by taking the increments in time and 
space small enough to guarantee that a sufficient number of stations exist for an accurate 
solution. A comparison test of the results for various increment sizes is an adequate 
means of determining appropriate increment sizes. 

SAMPLE PROBLEMS 
Spherical Cap With Dynamic Loading 

The first problem to be solved is one considered in references 12 and 13. An iso- 
tropic shallow spherical cap with clamped edges is subjected to a step pulse compressive 
pressure applied at t = 0 and removed at r = 5. The compressive nondimensional 
lateral pressure is taken as 60 percent of the classical buckling pressure p c ^ applied 
to a complete spherical shell where 



R 

and r n = — (fig. 4). The remaining shell properties are 

R = 100 in. (2.54 m) 

H = 1 in. (0.0254 m) 
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Ej = E 2 = 10 x 10 6 psi (68.95 GN/m2) 

<f> = 15.8094° 
v 12 = v 2 i = 0.3 

These parameters correspond to a shallow-shell parameter of A s = 5 where 
r -.1/4.V1/2 

X S = 2[3(l - |2Lj (31) 

and H is the maximum shell rise. The reference length CHAR or a is set at 100 in. 
(2.54 m) and E 0 is taken as 10 x 10 6 psi (68.95 GN/m2) with Hq and a set at unity. 
In addition, spatial increments are set at 1/25 of the meridian and the time increment 
(DTAU) is set at 0.25. The number of spatial stations and size of the time increment for 
this problem were established by comparing increasingly small spacings until stable 
solutions were obtained. The time response is desired out to r = 10 and a printout is 
requested at every fourth time increment. A complete listing of the program along with 
the results for this sample problem is contained in appendixes F and G. 

The namelist GIVEN quantities become 

N = 26 
NTYPE = 3 
CHAR =100. 

H0 = 1. 

S0 = 0. 

R0 = 100. 

PHI0 = 15.8094 
PP = -0.6 
PPS = 0. 

E0 = 10. x 106 
El = 10. x 106 
E2 = 10. x 10 6 
NU12 = 0.3 
NU21 = 0.3 
SIG0 = 1.0 
N0NLIN = 1 
C0NV = 0.001 
IDYM = 1 
KMAX = 40 
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DTAU = 0.25 
ALFA1 = 0. 

ALFA2 = 0. 

T1 = 0. 

T2 = 0. 

ITEMP = 0 
LBCL = 1 
LBCR = 3 
SL = 0., 0., 0. 

SR = 0., 0., 0. 

IFREQ = 4 
ISTART = 0 

The descriptive card comment is 

SAMPLE PROBLEM FOR A CLAMPED SPHERICAL CAP WITH LAMBDA = 5. 

The namelists ELB0L and ELB0R are not needed since neither LBCL or LBCR is 
set equal to five. 

The FUNCTION IP sets the time steps m at which there are step load changes. 
The input data (DTAU) was selected to have an incremental size, At = 0.25. Therefore, 
the abrupt or sudden load changes occur at stations m = 0 and m = 20 for r = 0 and 
t = 5, respectively. In the program the subscript m is represented by K0UNT. If 
there is no sudden change in load, IP is set equal to zero. If there is sudden change in 
load, IP is set equal to one at that time station. The user-supplied statements for 
FUNCTION IP(K0UNT) become 

IP = 0 

IF (K0UNT .EQ. 0) IP = 1 
IF (K0UNT .EQ. 20) IP = 1 

The user-supplied information to FUNCTION T(I) is 


T = 1.O/H0 


where T defines h at a station i-1/2. 

FUNCTION PL(I) and FUNCTION PS(I) set the lateral and meridional loads. For 
the sample problem a compressive uniform lateral load is used. Thus, in FUNCTION 
PS(I) the user-supplied statement is 


PS = 0. 
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Common statements provide the values for ETA ( 77 ), LAM (A), NU12 ( l '12)> 

NU21 ^2l)’ ROA (r Q ), PP, and K0UNT to be set in the function statement. Therefore, 
the user-supplied statements for FUNCTION PL(I) become 

PCL = 2 . *LAM*E TA/ (3 . * ( 1 - NU12*NU21))**.5*(T(I)/ROA)**2*E10 
PLL = PP*PCL 
PL = PLL 

IF (K0UNT .EQ. 0) PL = 0. 

IF (K0UNT .GT. 20) PL = 0. 

IF (I ACC .NE. 1) GO TO 1 
IF (K0UNT .EQ. 0) PL = PLL 
IF (K0UNT .EQ. 20) PL = 0. 

1 C0NTINUE 

Here IACC is computed internally at time points where the load changes abruptly (i.e., a 
step load) leaving the load doubly defined at that point in time. The first definition of a 
doubly defined load point at the mth time station (represented by K0UNT) is placed before 
the IACC statement card and the second definition of PL at that time-step point m is 
placed after the IACC statement. In this sample problem the abrupt load changes occur 
at K0UNT equal to 0 and 20. Therefore, at K0UNT equals 0 the statement before the 
IACC statement is 

IF (K0UNT .EQ. 0) PL = 0. 

and after the IACC statement is 

IF (K0UNT .EQ. 0) PL = PLL 

At time station m = 20 the statement before the IACC statement is 
IF (K0UNT .GT. 20) PL = 0. 
and after the IACC statement is 
IF (K0UNT .EQ. 20) PL = 0. 

At m = 20 the load is being suddenly removed so that initially the loading is equal to 
PLL and finally is equal to zero. The negative sign in the namelist GIVEN quantity PP 
makes the loading compressive. 

Since ISTART = 0, the FUNCTION DV(M,I) will not be called. It is interesting to 
note that since the initial conditions are zero, ISTART could be either one or zero. Since 
z i,0 = z i 0 = °> tIie in P ut to FUNCTION DV for ISTART = 1 would have been DV = 0. 

FUNCTION TDZ and FUNCTION TZDZ are completely contained as defined by equa- 
tions (16) and input constants for use in defining <*]_, Tj, T 2 , and j are provided 

for in the namelist GIVEN through ALFA1, ALFA2, Tl, T2, and ITEMP, respectively. 
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The statements appearing in appendix F for these subroutines assume that the reference 
surface is located at the mid surface. 

A listing of the program with this sample problem is contained in appendix F. A 
special nondimensional value A, the average inward deflection, is computed and printed 
for ease of comparison with references 12 and 13. The value A will be printed when- 
ever NTYPE = 3. The output to K0UNT = 4 is contained in appendix G. 

The results for 60 percent of classical buckling load (this sample problem) and 
other percentages are summarized in figure 5. The agreement with the results in ref- 
erence 12 is quite good. The discrepancies are attributed to differences in problem for- 
mulation and time increment sizes. The agreement with reference 13 is also good for 
P* = 0.4 but poor for P* = 0. 6 and r > 2. The discrepancy between the present 
results and those of reference 13 for P* = 0.6 is attributed to the use of a five-degree- 
of-freedom analysis in that study as compared with 26 finite-difference stations in the 
present study. Reference 13 reports a dynamic buckling load of P cr = 0.52. In refer- 
ence 12 the P cr is 0.65 which agrees closely with the present result of 0.68 as shown 
in figure 6. 

An extensive study utilizing the program for both static and dynamic buckling has 
been made in reference 14. Also contained in reference 14 is a thorough discussion of 
both static and dynamic buckling criteria. 

Thermally Loaded Clamped Cylinder 

This sample problem demonstrates the use of the program for analyzing thermal 
loads. The problem chosen is that one contained in reference 15 where an isotropic cyl- 
inder clamped at the first station and clamped and on rollers in the longitudinal direction 
at the final edge undergoes a linear temperature rise of 350° F (194.4 K) from one end of 
the shell to the other. The analysis is linear and the shell dimensions are 

Smax = 48 in. (1.2192 m) 

R = 12 in. (0.3048 m) 

H = 2 in. (0.0508 m) 

v 12 = V 21 = °- 3 

Ei = E2 = 28 X 10 6 psi (193 GN/m2) 

ai - = 9.5 X 10-6 in./in. /°F (17.1 X 10"6 m/m/K) 

Thus, namelist GIVEN becomes 

N = 21 

NTYPE = 1 

CHAR = 1. 

H0 = 1. 
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S0 = 48. 

R0 = 12. 

PHL0 = 90. 

PP = 0. 

PPS = 0. 

E0 = 1. 

El = 28. X 10 6 
E2 = 28. x 10 6 
NU12 = 0.3 
NU21 = 0.3 
SIG0 = 1.0 
N0NLIN = 0 
C0NV = .001 
IDYM = 0 
KM AX = 0 
DTAU = 0. 

ALFA1 = 9.5 X 10-6 
ALFA2 = 9.5 X 10~6 
T1 = 350. 

T2 = 0. 

ITEMP = 0 
LBCL = 3 
LBCR = 5 
SL = 0., 0., 0. 

SR = 0., 0., 0. 
IFREQ = 1 
1ST ART = 0 


The FUNCTION T(I) requires the statement 
T = 2./H0 

The FUNCTION TDZ(I) requires the statement 

TDZ = T(I)*T1*((FL0AT (I) - 1.5)/FL0AT (N-l))**2 

Since equation (17) is now simply T = Tj then Tj = 350^|^j 

The description card comment becomes 

THERMAL PROBLEM OF MENDELSON PG 186. 

Since LBCR = 5 then namelist ELB0R becomes 
OMEGAR (1,1) = 1.0, 8*0. 

ALAMDAR (1,1) = 4*0., 1., 3*0., 1. 

where the elements are read in columnwise. These input data correspond to the boundary 
conditions at i = n. 
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The right-hand side of this equation is supplied by the vector SR in the namelist GIVEN. 

Functions PL, PS, IP, DV, and TZDZ require only the following cards, respectively: 

PL = 0. 

PS = 0. 

IP = 0. 

DV = 0. 

TZDZ = 0. 

The results of this analysis are compared with those of reference 15 in figure 7. 

The comparison is good; the differences are attributed to the fact that in reference 15 the 
isotropic shell is approximated by a six-layered shell and the u deformations are 
neglected. 


CONCLUDING REMARKS 

A computer program has been developed to analyze thin shells of revolution which 
are both elastically and thermally orthotropic and are subjected to either mechanical or 
thermal loads. These loads can be applied either statically or dynamically. The program 
has many options concerning geometry, boundary conditions, and loading built into the 
subroutines. In addition, the basic subroutines of the program allow stiffness and geom- 
etry changes to be made easily without a detailed knowledge of the entire program. The 
present report describes the numerical analysis procedure and serves as the user's man- 
ual for the resulting computer program. 

A sample problem of the dynamic response of a spherical cap is included. The 
sample problem demonstrates the input data preparation for the program as well as the 
accuracy of the results obtained. A second example of a cylinder loaded thermally is 
included to show the input data required for that problem. Here again the agreement with 
existing results is good. 

Langley Research Center, 

National Aeronautics and Space Administration, 

Hampton, Va., September 18, 1970. 
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APPENDIX A 


REDUCTION TO SIX GOVERNING EQUATIONS 


The six governing equations comprising equations (18) are derived in this appendix. 
Equations (2) to (13) are reduced to six equations for the six unknowns njq, q, m-Q, u, 
w, and /3. Equations (10) and (12) are rewritten to define ejj and Kjj as 


D 11 

e ll = n ll 
A^G 


K 11 

G 


mu " N 12 e 22 " a - 2n 3 /< 22 + 1 1 

A^G 


K 


K 11 = 


11 


'll 


2 n ll + — m ll" 


M 12 e 22 


” Mo Key <y + • 


AG 


K 11 f m 
G 1 

(Al) 

K ll t n 

Ai 1 

(A2) 


When these expressions for e^ and are substituted into equations (11) and (13) 

the following equations result: 

n 22 = N 12 n ll + M 12 m ll + ( E 12 + C 22) e 22 + ( K 21 + K 22) K 22 + N 12 t l + M 12 t l^ “ ( A3 ) 


m 22 = N 3 n ll + M 3 m ll + ^ E 3 + ~^~j e 22 + ^ K 3 + “Jr)* 22 + ^1 + M 3 t ? 1 " 
where the coefficients are 
N 


(A4) 


11 - k 11 k 12) 


M 12=g(CllKi2-C 12 K n ) 

E 12 = 4;( 2C 12 k 11 k 12 - C 12 D 11 ' C ll K 12 

A Lt ' • 




K 21 =_ 2~( C 12 D 12 K 11 " C 12 D 11 K 12 + K 11 K 12 “ C 11 D 12 K 12) 

A. G 

«3 =4-(dhK, 2 - D !2Kn) 

m 3 = -|-(c 11 D 12 - K nK 12 ) 
a 2 g v ' 


(A5) 


(Equations continued on next page) 
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E 3=^f 


k 3 =^( 2Dl 2 K U Kl2 - D 11 k2 12 - C H D ? 2 ) 


(A5) 


G = J 2 ( C 11 D 11- k2 11 




The quantities e 22 and K 22 are eliminated from equations (Al) to (A4) by using equa- 
tions (7) and (9). Then substituting equations (A3) and (A4) into the equilibrium equa- 
tions (2), (3), and (4) yields the first three equations in equations (18). Substitution of 
equations (Al) and (A2) into equations (6) and (8) yields the fourth and sixth equations of 
equations (18). Finally, equation (5) can be utilized as the definition of /3 for the fifth 
equation of equations (18). Thus, the elements of the H, H, and M matrices and of 
the e vector in equations (18) are defined as follows. The elements of the H matrix 
are 


hll=^F^(l-Ni 2 ) 

h i2 = <p' 

h 13 = “ r m 12 

h!4 = -£2^±(e 12 + C 22 ) 

hl5 = - C ° S V‘ n0 ( E 12 + C 22 ) 

hi6 = co^^2i + K 22 ) 

*21 = -(<*>' +^N 12 ) 

h 22 = co = * 

. sin cp 

, cos cp sin <p/ 0 

h 23 = r M 12 

h 24 = - ^.2 ( E 12 + C 

. sin^rf)/ _ \ 

h 25 = ~ r 2 fl2 + c 22) 

cos <p sin cp/— 

*>26 - l 2 (K 2 l + K 2 

, cos <p XT 

h 31 = r N 3 

h 32 = ~ x " 2 

h 33=- E2 |^( 1 - M 3) 

h 3 4 = -2^4 (e 3 + x- 2 K 22 ) 

h 35 = - COS V‘ n *( E 3 + "- 2 k 22) 

h 36- £2 ^( K 3+ x ’ 2l) 22) 




(A6) 


(Equations continued on next page) 
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h __ D u 

41 X 2 G 

1143=-^ 


h 4 5 -*' + N 12 SBJfc 


h 51 = h 52 = h 53 = h 55 = 0 
h 56 = - 1 

h62 = 0 


, cos d> 

h 64 9^ m 12 

rX^ 


, COS 0 
h 66 = —7— m 3 


The elements of the H matrix are all zero except 

hll=f£ 


h 42 = 0 


h 44 = 


COS 0 -|^ T 

= r N 12 


u .9 COS 0 XT 

h 46 = N 3 


h 54 = 

h K H 

1,61 = -35 
C 11 


h 63 = 


h 65 = 


sin 0 


rX 


M12 


~ cos 0 a p' 

~ n ll 




h 26 = ' 


V 




The elements of the e vector are 
cos 0 


ei = 


N12*; + M U tf - tgi . Ps 


sin 0/, T ,n ,, ,m + n\ 
e 2 = — jr J£1 (Ni2t 1 + M 12 t 1 - t 2 J - p 

cos 0 /« T ,n ,m .m\ 
e 3 = — F^V N3t l + M3t l ' *2 ) 


(A6) 


(A7) 


(A8) 


(Equations continued on next page) 
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The only nonzero elements of the mass matrix M are 
m 14 = m 25 = 1 


(A8) 


(A9) 
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TIME DERIVATIVES AND INITIAL CONDITIONS 


It is seen from equations (18) and (A9) that the time derivative terms arise only in 
the displacement vector y^. Thus, as shown in reference 7, the y^ m terms can be 
written as 


y i,m 


1 h. 


(At) 


2\°fei y i,in + ^m^m-l + 



(Bl) 


where m indicates a time step. The constants a m , /3 m , y m , and "6 m are defined 
in reference 7 by use of Houbolt's initial starting procedure for homogeneous initial con- 
ditions of Yi o = h o = °- The P rocedure is presented here for general nonhomogeneous 
initial conditions where yj q and y^ 0 are given. From equations (2) and (3), q 
can be calculated. Also the following difference equations at m = 0 can be used to define 
Yi,0 and Yi,0 as 


%0 



1 

6(At) 

1 

(At) 2 


[2 y i,l + 3y i,0 * 6y i,-l 
( y i,l - 2y i.O + \-l) 



(B2) 

(B3) 


These two equations can be rewritten to define the fictional time points y. j and y. ^ 
Then use can be made of the general backward-difference equation 


y i,m 


1 

(At) 2 



^ y i,m-l + 4y i,m-2 



(B4) 


Therefore, by using equations (B2) and (B3) to eliminate the fictional points from 
equation (B4) at m = 1 and m = 2, values of c^, 0 m , y m , and 6 m of equation (Bl) 
are obtained as follows: 
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At m = 1, equation (Bl) becomes 







(B5) 


and 


Q!i = 6 /3j = -6 = 6j = 0 

At m = 2, equation (Bl) becomes 

" i ’ 2= ^( 2i ' i ’ 2 ' 4Si - 1 + 2yi '°) 

and 

012 = 2 ^2 = Y 2 = 2 

At m = 3, equations (Bl) and (B4) are identical and 

% = ^ ^m = = ^ ^m = 

Equations (B4) to (B6) completely specify all time derivatives and initial conditions. Sim- 
ilar results for nonhomogeneous initial conditions were obtained in reference 16. 

This procedure of using initial conditions is employed at every time point where 
there is a sudden change in load. 


- y M (B6) 

62 = 0 
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DEFINITION OF BOUNDARY CONDITIONS 

The system of equations defined in equations (18) and appendix A requires three 
boundary conditions at each edge. These conditions are derived in reference 4 and are 
defined by a combination of the following variables: 


r 

nn 


fu l 

< q 
mn 

> or < 

w ' 
/3 



J 


Thus, z is a six-element vector and the boundary matrix contains only three equations 
at each end. Therefore is divided into two subvectors x^ and y^ where 



as shown in figure 3 the left boundary is at station i = 1 and the right boundary is at 
i = n. Therefore 


DjSxj + Aj6yj = Zj 


°n 6x n + A n 6 y n = *n 


(Cl) 


/ SrnojA 

where the subscripts 1 and n refer to the first (s = 0) and last Is = — - — -I stations, 

respectively. In the program the vector l at i = 1 is SL and at i = n is SR. Both 
SL and SR are three-element arrays. The following conditions can be applied to either 
boundary: 


For a pole point where u = q = j3 = 0, the matrices 


S2, A, and l become 



— 


— 


— 


— 


r 


0 

0 

0 


1 

0 

0 


0 

n = 

0 

1 

0 

A = 

0 

0 

0 

l = < 

0 


0 

0 

0 


0 

0 

1 


0 





_ 


— 


— 


L J 


(C2) 
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For a pinned boundary where u = w = mu = 0, the matrices O, A, and l become 


D = 


0 0 0 
0 0 0 
0 0 1 


A = 


10 0 
0 10 
0 0 0 




l =< 0 


0 


(C3) 


For a clamped boundary where u = w = /3 = 0, the matrices O, A, and l become 


O = Null matrix 


A = I 


l = Null vector 


(C4) 


For a free edge where njj = q = mjj = 0, the matrices become 

fi = I A = Null matrix l = Null vector (C5) 

Finally, for a boundary with general elastic constraints, D and A must be 
defined by the particular problem and read in through namelists ELB0L and ELB0R. 

The vector l is always read in through namelist GIVEN. For these boundary condi- 
tions, the elastic boundary conditions on the left (i = 1) edge are read in through namelist 
ELB0L and on the right (i = n) edge through ELB0R. In other words, all nine elements 
of both S2 and A must be specified. 
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DEFINITIONS OF A, B, C, D, E, AND q MATRICES 


Equations (22), (26), (27), and (29) are related in the following manner: 


f k - F i~ 1/2 


+ G i-l/2< yi > " L i-l/2 


where k = 1, 2, . . ., 6 and i = 2, 3, . . n. For k = 1, 2, 3 


A2 (i - 1) = Nu ^ matrix 


2 (i - 1) ~ i = l F i-l/2 + 6 2,j + l 


(j = 1, 2, 3) 


. 8f k _L +6 ( n ll^)j-i/2 [ n l 1 /c°s 0 1^ 

'2(i-l) i-!/2 1.5-5 2 ~r, °2,j-4 V \ 2 v As 


0 = 4, 5, 6) 


D . JJ±J G . 4. s . w 

2(1—1) 8xj y 1— 1/2 2,j+l f 


(i = 1, 2, 3) 


F _ = J G . 6 ^ n ^i-l/2 + 6 pi 1 /cos + _1_' 

E 2(i-1) s yi S G i - 1/2 l,j-5 2 r\ 2,j-4 7) ^ 2r As 


• 1/2 Ai 


^2(1-1) = Api-i-yi-i’^i’ 8 ) 


(j = 4, 5, 6) 


where the barred vectors indicate the approximate solutions. For k = 4, 5, 6 




B - 8f k _/ F . 6 . hll/1 

21“ 1 8y^ j y i-1/2 4 , 3-2 4^ ^ 


(j = 1, 2, 3) 


(j = 4, 5, 6) 


(Equations continued on next page) 
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af 


C 2i“l = a Xi " ( G i-l/2) kj 


r, _ 8f k _ ( n * Pi- 1/2 

21-1 -Jf ~ y i-1/2 - 6 4,j-2 


^2i-l = Null matrix 


q 2 i-i = - ^(*i-i»y i . 1 » x i»y i » 8 ' 


kj 


"N 

(j = 1, 2, 3) 


(j = 4, 5, 6) 




(D3) 


J 


Equations (Cl) are related to equations (24) and (29) in the following manner. For 

i = 1 


Ci - 
Di = A l 

> 

Ej = Null matrix 


(D4) 


11 =*L 


For i = n 


A 2 n = Null matrix 


®2n - 


C 2n " A R 


q 2n = l R J 


(D5) 
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The complete set of equations becomes 


“ C 1 

D 1 

E 1 

F 1 




X 1 


r *1 " 

B 2 

c 2 

»2 

e 2 




yi 


q 2 

A 3 

B 3 

c 3 

d 3 

e 3 

f 3 


x 2 


q 3 


A 2(i-1) 

B 2 (i— 1) 

C 2 (i — 1) 

D 2(i— 1) 

E 2(i-1) 

< 

*1 


q 2(i-l) > 



A 2i-1 

B 2i-1 

C 2i-1 

°2i-l 

E 2i— 1 

yi 


q 2i-l 




A 2n-1 

B 2n-1 

c 2n-l 

°2n-l 

x n 


q 2n-l 

— 




A 2n 

B 2n 

c 2n 

yn 
V n J 

\ 

q 2n^ 


where Fj and Fg are null matrices added for a later programing convenience as 
discussed in appendix E. Since they are null matrices at this point, they do not affect 
any of the previous definitions in equation (29). 
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RECURRENCE EQUATIONS 

If for convenience 6xj and 6yj are represented by Xj and yj, a recurrence 
solution to equation (D6) can be obtained based on the Potters method (ref. 9). To insure 
nonsingularity, elements Cjj, C22, and C33 of matrix Cj must not be zero. If 
either cjj or C33 is zero in Cj, then row one or three of Cj_, Dj, Ej, Fj, and 
qj must be interchanged with row one or three, respectively, of B2, C2, D2, E2, and 
q2» If C22 is zero, then several row manipulations must take place. First, row two of 
Cj, Dj, Ei, Fi, and qi is placed in row two of B2, C2, D2, E2, and q2; row two 
of B2, C2, D2, E2, and q2 is placed in row three of A4, B4, C4, D4, E4, and 
q4; row three of A4, B4, C4, D4, E4, and q4 is placed in row three of A3, B3, 
C3, D3, E3, and F3; and row three of A3, B3, C3, D3, E3, and q3 is placed in 
row two of Bi, Ci, Di, Ei, Fi, and qi to complete the cycle. Simpler substitu- 
tions could be made for specific shells and boundary conditions but the preceding row 
interchanges yield nonsingularity for all shells. Elements can exist in the Fi and F3 
matrices of equation (D6) as a result of the row interchanging. Since C 1 is nonsingular 
then 


x i = p iyi + Qi x 2 - c I lF iy2 + R 1 


(El) 


where 


Pl = -C^Dj 

Qj.-q'Eil 

R i = c iV 

y 

and, in general 


x i " P 2i-l y i + 92i- 1*1+1 + R 2i-1 
yi = P 2i x i+l + Q2iy i+ 1 + R 2i 


(E2) 


(E3) 
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where 


G r - |^A r P r _2 + B r jP r _i + A rQr-2 + G rj 
P r = G r J^A r P r _2 + B r ^Q r _j + D Jj 




Q r _ _ G r E r 


R r = Gj. 


^r “ 


r tt r-2 " ( A r p r-2 + B r y 


AvR 


:) R r-l] 


The exceptions (in addition to eqs. (El) and (E2)) are as follows: 
At r = 2 


Q 2 = -G; 1 ^ - BaC^F^ 

At r = 3 

P 3 = P r h-G^AsC^F! 

x 2 = p 3 y2 + Q 3 x 3 ■ G 3 lF 3 y 3 + r 3 

At r = 4 




Q4 = Q r + 

At r = 5 


(a 4 p 2 + E^Gj'Fj] 


p 5 = P r + G^AsG^Fj 


(E4) 


(E5) 


(E6) 


(E7) 


(E8) 


Also at r = 2N, Aa^, DaN> E 2N» an( * e 2N-1 are null matrices and equations (E3) 
reduce to 


y n = R 2n 




x n = p 2n- l y n + R 2n-lJ 


(E9) 


Thus, the procedure of solving equation (D6) is to sweep down the diagonal to solve for 
y n and x n and then back up the diagonal using equation (E3) to solve for all Xj and y^. 
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The 3x3 matrix multiplications, additions, and inversions are performed by the CDC 
library subroutine MATRIX as detailed in reference 17. 
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PROGRAM LISTING 

The program listing for the sample problem is as follows: 



P H* 0 ' ■« 

v \ - 


^ A ) A U N ( i • K 0 r • 

03 i PU \ , T APtD= INPUT ) 


c 

iuy 4= i > o--i i r 

"} 

1 “i ^ r ■< m hi -> j t. N t 

KHnPunSF. Pk OELO i 


c 

lur^=i INC. 

01)** 3 T-cl T k*- 1 

c. >J 1 PCbK^J^4^t. PKu 3 Lc^ 


c 

I 0 Y Li = 

T J E 


STATIC r, JC^ L I 

i- LO^U ''.ILL 3 L CaLCUL^iEl) 


c 

IN 

T ■ i 


V? A » I u '■'i -MU'* 0 r. ^ 

ij h I Ii'-iL sTcPS WrlE.M loY i ^ = 

1 • 

c 

< •■* A X I - 

T >i: 



*j r *5 T a r I C c 5 ULNLl rv ib STEPS 

WHC.N 1 UYN 

c 

- T Y P F •- 

0iJ»' 

L 

r •• x i * cyl 

1 1 ,UC K 


c 

T Y p F t- 

0‘ 1 ‘ 


T ./ " j N - '-I; \ 

11 


c 

• T YP p F 

\»0 

l_ 

T 1 j IN 1 ^ u ^ 

H > v X 


c 

>rypp F 

J'j 1 


T ) -t- 1 Hl OH 0 ’'* 

?. 1 k y ,rf lUbT bPtClFltJ) IN 

Gto-irr. 

c 

I p L-'CL 

= 

X 

T -r. L c F i rt 1 

M J vj L r pul N T 


c 

IP L <C.L 

r: 


T -1 c_ L c F i " C i 

u i vjiMLU 


c 

IF l^CL 

“ 

;> 

f He LcF 1 r>C i 

''hi ^cu 


c 

if l * Cl. 

- 


T LtiF • ' N T 

- F-hL 


c 

IF L H C L 

= 


1 i : ! L-r • -'C 1 

^ a 1 1 a 1 - 1 i oF llaSTIC CONTACTS 

c 

IF L -CP 

- 

X 

T F r* C 

1 " -- p 3 L t po 1 i'j r 


c 

IF L.L-CP 

- 


Ti- -JbM -C 

L “•» K I 


c 

IF L-C- 

r 

; 

T 1 1" w i 0 n I -C 

IS r 1 ALU 


c 

IF l. *-< C p 



T’ll ^ l'.;nl -«C 

1 - r u*.t 


c 

IF 1 . k 0 ^ 

- 


T - J 3 . -3 » 1 C 

i ^ >3 ■ J *■< i «: I X OF t l A b T 1 C C ON S T a j T S 

c 

1 F MO-i. 

I - 

- 

•a L i C ** ■ n e)I 

•J f US i ^ 0 ri 1 A i ,\J c U 


c 

IP is-O-t, 

i *■' 

- 

; * v‘i vL 1 •- -p 

s u l 1 j i 1 0 v j 13 U t 5 T A 1 N E U 



IF' IM'-r r ‘ I -'-'L.- 'ILL nt -(DLVLU FUR I HE STATIC L0«us FIRc,T 

IF J -> T •> ' f =. 1 '111 H >1 1 I-. -ita A"gu VELOCITIES Mu.bT HE uIVFN 

C O v * ' j , / - l 1 / *0 * * • • k n f :> u a • Ob* Ok 

(. > ) • A / .i L /.«/ 'i* *:>-l (-)•>)■ h l A|.-i jAL(iiJ) » 0 *4 L b A ^ ( j • J ) , ALA M DAP (3*3), 

i - 0 ( 3 ). .* ( * J 

cn v.F'j ./r.^ /u.. , i Cv., I .11 

Or •'!/./ l 

C o ’ - s > 1 / L 3 /' ( - v ;« j» j) , ■ • (._ 3d • i) » l < 3,3) » 1 ; { u d » J t 3 J *E(E0E,3,3) , 
1 > ( Z M 

C o i • *■ / ■ l 1 / C 1 i * C ± ' • viic., r\ x i ♦ id * ft E E , D 1 i * L>1 E , Odd 

C ) v * v > 1 / ■ L ( / .* i |_ • j ( -j 1 1_ , 3 ) • c- k 3 c 9 t p b 

Co v‘ [.,(*■ v, / *. L / L ~ * 1 « ri‘t -c‘** oU^if t TM 

C < ') ■ * * 1 / : - L • / I * j t ‘*i , * 1 * « • 4 u \ L i 4 , L ; 1 1 
C 'i‘j j / ", j 0 / ( i • ' I ) m -* 1 l x i ) * 0 P ri i ( i 0 1 ) 

C'li * ' J %' / ■' \ 1 1 / ‘A F* a x * a l r ^ £ • lx, T d * 1 1 t V 1 P 

C <V : ‘‘ ’’ / L t : / ' 'j ) 4 1 

C 0 v J • 4 / i_ ; )/ A i(i •'!*-)■ % L ( A O 1 « ) , X 3 ( 1 J 1 • d ) 

0:) v ’ f 1 / .il. t 4/ b < • L « ■* ►" “ J 1 , *a \J 

CO 'i ■'•;»••• / Li >/ * ; 

C 0 v ' ' / ' N / ' ’ 1 _ i . / i p K * 1 
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D I Mt NS I UN XI2C2.3) 

DIMENSION X B ( 202 » 3 ) 

REAL LAMS 

REAL MJl 2 * N021, Kll, K12, K22 

REAL N12, F3, M12, M3, KR12, K3 

REAL LAM 

1 FUR MAT (12 X, *TFE LEFT OMEGA MATRI X*, 40X*THE LEFT LAMBDA MATRIX*, 

1 //5A*(* 3IE12.5, 5X ) » * ) ( N ) +(*3(E12.5, 5X)*MU) = * E12.5/ 

2 5 X * ( * 3 ( E 12 • 5 , 5X ) » * ) ( Q ) +(*3(EI2.5, 5X)*)(W) = * E12.5/ 

2 5X* ( * 3IE12.5, 5X),*)(M) +(*3(E12.5, 5X)*)(B) = * EI2.5) 

2 FURMAT ( 12X, * T FE RIGHT OMEGA MATR I X*, 40X*THE RIGHT LAMBDA MATRIX*, 

1 / /5 X* ( * 3 ( E 1 2 . 5 , 5X»,*MN) +<*3<E12.5, 5X)*MU) =* EI2.5/ 

2 5 X * ( * 3IE12.6, 5X»,*)(0) +(*3(E12.5, 5 X ) * ) ( W ) =* EL2.5/ 

2 5 X * ( * 3IF12.5, 5X ) , * ) ( M ) +{*3(E12.5, 5X)*)(B) =* E12.5) 

b FURMAT (10X , 15, 2(4X, Ell. 5), 2(9X, 16), E15.5) 

f FORMAT (EOF 

1 ) 

13 FORMAT 1/5 X ,*S TAT ION NO*5X*U D0T*15X*W DUT*15X*U DUTDOT*12X*h DOTDO 
IT*/ / ( 5 X , I 5, A ( 5X * E15.8) )) 

30 FURMAT ( 3U X * I F IDYM=0 THF SHELL IS LOADED ST A T I C AL L Y* / 3 0 X* I F IDYM=i 
IT HE SHELL IS LEADED 0 Y N AM I C AL LY */ 3 OX* I E I0YK=2 THE SHELL STATIC BU 
2CKLING LOAD IS C AL CUL AT ED* // 30X *F JR THIS RUN, I DYM= * I 2 0// 3 OX* NUMB E 
JR UE STATIONS =* I 1 9 // 3 OX*MER I C I AN/ REF ER ENC E RATIO =* F20.6 // 3CX 
A * R A D 1 US/ REFERENCE RATIO =* , F2 2 . 6/ / 30X , * THI C K/ REF RAD RATIO =* 

3 F2A.b//15 X*T 1 / F 0 =* F20.6,23X*E2/F0 =* F 2 3. 6// 1 3X *NU1 2 =*F21-6 

o,2 3 X*NU 2 1 =*F24.6//13X*CU/S0 =*F20. 6 , 23X*REF DIST =*E 2 0 . 6/ / AOX* I F 

5NUNLIN = 0 CNLY L INEAR TERMS ARE. US ED* / AO X* I F NUNLIN = I NUNLINEAR 
6 TERMS US t' G * / / A 0 X *F G T THIS RUN NUNLIN = *16/) 

03 F URMAT {//AX, ♦STATION NU*,7X,*N-S RESULTANT*, 

1 5X,*SHEAR FORCE* ,7 X , * M- S RESULTANT*, 5X, *U-D FT URM A T ION * , 5X, 
2*W-DtF0u,.|AT ION*, SX, *Bt T A ROTA T I UN*/ / ( AX , I 6 , 6X , 6 ( 3 X , h 1 3. 7 ) ) / ) 

9/ FORMAT (lLX*M:. OF ITERATIONS= *16, 19X*ERR0R NORM =*E15.8, 1CX*CYC 
1 L E 3 THIS ITERATION =* I 6 ) 

203 F URMA T ( hX,*PL(I) =*C I 5 .3, 3X , *PS( I) =*E 13 .8 , 3X*SL ( 1 ) =*F15.fl,3X 

1 * 3 L ( 2 ) =*E 13.8, 3X*SL ( 3) =* E 1 5 . 8/5 X* LUA C CHANGES =*Ib,5X*DEL =* 
2F9.6.S X* I Tf RAT IONS =*I6,AX*N0 UF CYCLFS = * I A , A X *XNOKM =*E12.5) 

117 F ORMAT (AoX, * T F E MAXIMUM SHE1L RISE IS *F17.3) 
llo t- URM AT ( A o X , * T F C SHELL THICKNESS =* T17.3) 

231 FORMAT (SX*Rf'AX= *F15.3, 5X*I= *16) 

z-v i format <//acx*tee geumftky uf the shell folluws*/5x*station ng* 

1 AX* MERIDIAN LIS* 8 X*R A E I AL D I S *1 0 X * ANG L E ( R AD ) * 10X* CUR V AT UP E ( CF I / D 
2S )*/ ) 

2 A2 FORMAT (‘3X, 16, A(5X, E16.6)) 
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t- LRiXAT (/I CX*KCLNT* 10 X* TAli* 12X*DCL* 12X 

1 YCL H S=5 AX *XNURM*) 

24-*+ Fu^ATI-l*) 

.i. Ai v i tL I S T / C IVFN/NTYPF , N » Rf* SO, HC, fcO, 
li'iDXI li\, CCNV, PRIC, LiKL, L B C R , PP, PPS, 

2 IOYM, F RAX, C TAtJ , ri, T2, ALTAI, ALFA2 
iMA 1t.LI SI/ El’M'L/ UMtGAL, ALAMCAL 
NAMELIST/ ELBCP/ (j ME GAR , ALAMOAR 

DATA PI / 15^2 t>b 3 b 8979 i / 

1 u 1 k t A U b I V h N 

IE ( L-CF, SI 3 C l , 302 

301 SHP 
TJ2 UiMTINUE 

PRINT GIVEN 
READ 7 
PRIM 7 

IF (LUCL .ft. 5) REAL) CLRTJL 

IF (LHC' . fcC. 5) EFAU EL3GK 

LA'-1= HC/LHAK 

3UA = Sl')/(.F> AP 

,< L A = R U / C r> A R 

LIAMU/31GC 

tlO= tl/EL 

fc 2 0 = t 2 / fc L 

Ri; 1C.'LVV=1 

IbUCK.= 2C 

K. S T U P = 0 

I ACC=0 

0 A V fc = 0 . 

1 N I T = 0 
KuUNT = 0 
ITFk =C 
LLH ANGE=C 
XNUPM=0. 

DP = PP 
L)P >=PPS 
OS l = S L ( L ) 

DS22=SL(2 ) 
u i> 3 - i> L ( 3 ) 

0 T 1 = 0 . 

L)T 2-0 • 

PRINT 2 44 

CALL GfcCPiT V (NTYPE ) 


ITERATIONS* bt 


I, fc?, NU 12, NU2 l , SIGO, 
SL, SR, CHAR, 

ITtMP, IFREQ, ISTART 
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PRINT 6C, IOYM, N, SNA , RG A , LAM, fc 1 0 » E20, N!J12, NU21, ETA, CHAR, 
1 NO.NL I N 
US2=DS/2. 

SMER=— CS 

DEL = ( PH I (2 I- PH l C D ) /?. 

PRINT 241 

DU 240 1=1, N 

RI= k { I )+DS2*CCS(PHI{ I ) ) 

FI = PH I (I J+CEL 
b«cR= S ME K + C S 

24U PRINT 242,1, S^FP , RI, FI, DPHI(I) 

IFINlYPt. NF. 1 ) on TfJ 211 
2U = R CA-t- ( 1 .-CP'S ( RPHI ) ) *FTA 
2 SUP = /u*2C*REA/ (ETA* 2.) 

Zti A= ZC/EIfi 
PRINT 117, ZCA 
211 CIJN TIN Lb 
I T = T( 1 ) 

PRIM lib, TT 
L = 2*N 
L 1 =L- 1 
C PS = DT AU» 

FPS2=OTAU**2 

Du d 1 J = 1 , 1 
D C dl M = 1 ,? 

XULCH J , M) = C . 

X B f J , M ) =0 . 

2 1 X ( J , Fl ) = t . 

DO 2 c 3 T = 1 , N 
Ui: 2b M= 1 , 2 
X 1 ( I ,.A )=0 . 

X 2 ( I , M 1 = C . 

2 J A J ( I , M ) = u . 

I TEK=0 
NCYCLfi = 0 
TAU=-r- PS 
id HE! = I T L P + 1 

hC YC LF = Nu YC L F + 1 

CALL RCCuN (LHCL, LBCR, X, ITER) 

H (ITER . NF. . 1 ) Gf! TJ 11 

PRINT 1, ( (C< l,M,MM) , MM=1,3), (0(1, M, MM, MM = 1,3), SLIM), M=l,3) 
PRINT 2, ( ( i<( L, M,MM! ), ^,= 1, 3), (C(L,M,MM), MM=l,T), SR(M), M=l,3) 
11 CONTIlUf 
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IF ( I STAR I .EQ. 1 J CALL INITIAL (X, LBCL) 

IF ( ISTART .Nfc. I I GO TO 10 
IF ( I TER .EG. U GO TO 93 
10 00 84 1= 2 »N 

CALI A BtDE S ( X , I, NCYCLEJ 
84 CONTINUE 

CALL FCNCTIX, N) 

CALL PCTTER(X) 

IF ( ITEK .EC. 1) GO TO 94 
XNOPM^O. 

DXSUM=0. 

00 9b J=1,L 
DU 95 y= 1,3 

XNOKM= X NOR N+X CLO { J , P) **2 
95 DXSUM= DXSUN+X ( J » M ) * $2 

IF (XNOKM .EG. 0.) GO TC 94 
xNO RM= SORT (CXSUM/XNORM) 

9 h 00 9b J=1 ,L 
DO 96 P = l,3 
9b X ( J , M ) =XOLC(J,P) 

IF (ICYP .EG. C) 

1PKINT 85, (1, (XI 2*1-1 , J) , J = l,3 ), ( XI2*I ,J ) , J=i,3) , I=1,N) 


i r 

( { I CVM 

.Ei J . 

0) 

.AND. ( iNiCNL IN . EQ. 0)) 01 TO 105 

if- 

( ( 1 tYfV 

.f 0 . 

0 1 

.AND. f ITEK .EQ. II) GO TO 92 

IE 

( XNuKfo 

. EQ . 

0. 

) GU TO 9 3 

IF 

( I T tP. 

.rc. 

1 ) 

GO TC 92 


IF ((NCYCLfc .GT. I8UCK) .AND. ( IOYM .EC. 2 )) GO TO 201 
IF ( XNCRN-CCNV ) 93, 93, 92 
9 5 CON TINLE 

TAO=T AU+EFS 

IF ( MCD( KCU.'T, I FREQ) .Nfc. 0 ) GO TO 9 

8 PRINT 85, II, (X( 2*1-1, J) , J=l, 3), ( X( 2*1 ,J ) , J=l, 3) ,1=1 ,NJ 
105 CALL ANSWLRSINTYPE, X, L8CL, LBCR ) 

IF (IP(KCUNT) -EO. 1) GOTO 103 

IF ((ICY! .FG. 1) . AND. ( MCC ( KOUNT , I FRED ) .EQ. 0)) PRINT 243 

105 CONTINUE 

IF ( 1DYP .EG. 0 > GO TO 101 

9 CONTINUE 

l F ( NTY PE . N E .3 ) GO TC 206 
wSUM= C. 

DC i 5 1 = 2, N 
J= 2 * I 

5 WSUF=wSUfc+(X( J,2 ) +X (J-2,2) ) / 2 . *R ( I )*CUS(PHl ( I))*DS 
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20 o 


4 

12 


2 4 
2 o 


DAV E = W SUM / Z SUN 
CONTINUE 

I F ( IDYN.tC-2) GO TO 204 
IF (ISTART .EC* 1) GO TO 12 
IF (1PIKGUM) *NE ♦ 1) GO TO 3 

DO 4 J = 4 , L , 2 
DO 4 M= 1 1 2 
I - J /2 

XJ=< XI J,K)+X( J-? ,M) >/2. 

XD =< ll.*XJ-L8«*XL( l v M)+9. *X2U f M) -2. *X3< I , M ) ) / ( EPS*6 . ) 

XI ( I » M ) =X J 
X 2 ( I , M ) = XD 
I ST ART = 0 
NC YCLE = C 
INI T=KOUNT 
I ACC = 1 
DU 26 1=2 

CALL ABCDFS (X * I, NCYCLF) 

J=2*I- 1 
DO 26 1 » 2 

A X= 0 • 

BX = 0 * 

cx=o. 

DX = 0 • 

ex=o. 


DO 24 Rf'=l 
AX = AX+A(J 
bX= rtX+EHJ 
CX= CX+C(J 
L) X- D X + L ( 2 
IF (I . LG . 
t X= LX + h ( J 
COMT INUb 


, * ) *X( J-l , Mi*) 

, )*X< J ,MMJ 

♦ *■ )*X< J + l t MM ) 

N ) CO TO 24 

♦ F )*X ( J + 2 , MM ) 


XJUt*)=AX*FX4CX + CX*EX-S(JfM) 

PRINT 14 , U, X 2 ( I » U » X 2 ( I » 2 ) » X 3 ( T t l ) t X 3 (I, 2 ), 
PRINT 2 a 4 
1 A C C = 0 
COMT l N U F 

IF { 1 D Y i v .bC* 1 ) PRINT 6 f KCUNT, TAU, DAVE, ITtK, 
tviCYCLF =0 
KUUM=KL CRT + 1 

IF (KiJUNT .IQ. (KMAX+ 1 )) GO TO 102 
I F- ( IiJYF - F L . 1 > GO TO 92 


1=1 tN) 


NC Y CL E, 


XNORM 
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2 3<t CCijriNUF 

if- (IUY* • N F . 2 ) GC TO 102 
Pi: INI Ml, 11, T 2 

Vi i FORMAT {Va, *Tl = * F I V • 7 , 5X*T2=*F IV. ?) 

WMA A = 0 • 

01. 2 '30 l-l fN 
J = 2*I 

if ( X ( J 1 2 ) ME. W MAX) Gu TU 230 
I C •— I 

WMAX = X { J T 2 ) 

20u CUNTIMt 

PH\T 22-1, MAX, IQ 
MUM =KLlf\T + l 

IF ( ,\i I Y P r_ .IM. 3) DAVE- NMA X 

Pf.IM 2u* , FF, PPSt (SLID, 1-1,3), LCHANGC, UAViz, ITER, 
i XoU R M 

IF ((NTYPt •FC. 3) -AND. (ABS(DAVt) .GT. 2-0)) GO TO 101 
20 i DNE=1. 

IF (K-;)U\T -GT. KM AX ) GO TO 102 

I F (LOANGb . f C. 3 ) 0 N F =— 1 . 

pp = PP + l)P*CN t 

f PS=PPS+JFS*CKt 

S L ( 1 ) = SL( 1 ) *0$ 1-fJNF 

SL( 2 ) = S L ( 2 )+CS2?*CNE 

SL( 3 ) - S L ( 3 )+nS3*GKE 

) 1-n + cn *cm 

T2= T2+012 MKF 

IF ({kSTuP .hC. ( KHUN T — 2 ) ) .AND* (LCHANGt .CO. 5)) UNE=-. 
IF ((KSIuP .re. { KGU N f-2 ) ) -AND. (LCHANGE MO- 5>) GU T(J 
NUN CON V = 1 

if INCYCLF .LE. I RUCK ) GO TC 202 

DCi 20 V I-1,N 

J = 2*I— 1 

DU 2 0V iS - 1 i 7 

XCL 0 { J , 1 ) - X 1 ( I , M ) 

XUL D { J ,2) = X2< I ,P) 

XULIM J,3)=X3( I,M) 

2 0 V J-J + l 
22o C. 0 N T I N U C 

KSTOP= KL’IM 
NUN CUN V = 0 
PP=PP-CP* 2. -CNF 
PPS=PPS-UPS*2 -*CNE 


NCYCLF, 


5 

226 
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SL11 )=SL( 1 )-0Sl *2.*0NE 
SL( ?) = Si ( 2>-C$?2*2.*QNE 
SL ( i) = SL< 3 ) -OS 3 *2.*ONC 
T i = Tl — CTl * 2 * * C N E 
^ 2=I2-CT2*2.*CNE 
OPS=DPS/5 . 

DP= OP/ 5. 

US 1 — DS 1 / 5 . 

D S2-D S 2 2 / 3 . 

US3 = DS 3/ 5 • 

DTl=DTl/6. 

U T 2 = DT 2 / 5 . 

LCHAIMGE = LCMNGF*1 
If {LCHANGE.EQ.5) GO TO 215 
if (LCHANGE .EG. 6) GO TO 101 
GO TO 224 
215 CONTINUE 
PP= PP- D P 
PP$ = PPS — DPS 
SL ( 1 ) = S L ( 1 )-CSl 
SL ( 2 ) = SL ( 2 ) -0 S 22 
SL ( 3 ) = S L ( 3 )-CS3 
T L=ri-UT1 
1 2 = T 2 — D T2 
2Ud DO 216 1=1, N 
J=2*I-L 
DO 216 M — 1,2 

X ( J , 1 ) = X1 ( 1,M)+5.*(X1( It'D — XB(Jfl) ) 
X( J ,2) =X2 ( 1 ) +5.* (X2 ( I J ,2 ) ) 
X ( J t 3 ) = X 5 ( I , V ) + 5 . * ( X3 ( I , > - X B (J , 3 )> 
XC L D { J , 1 ) =X ( J , 1 ) 

XUL D { J , ?) = X ( J ,2) 

X ( J L D ( J , 3 ) - X { J , 3 ) 

2 l o J=J+1 
224 CGNTINLF 
PP = PP + D P 
PPS-PPS+UPS 
SL ( 1 ) = S L ( l)+nsi 
SL( 2) = SL12 MCS22 
S L ( 3)=SL(3)+CS3 
I 1 = T 1 + C T 1 
T2=T?+CT2 
202 NCYCLE = 0 
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DO ?07 1=1, N 

J=2*I— 1 
DO 207 M= 1 , 2 

IF (NONCONV .EQ. 0 ) GO TO 209 

xb( j, i >=xi ( i 

X B ( d , 2 )=X 2 ( I, K ) 

XB ( J,3)=X3(I,P) 

X 1 ( I , M ) =X ( J , 1 ) 

X2 (I ,M)=X ( J,2 ) 

X3( I , M J=X ( J , 3 ) 

209 X ( J , 1 )=XLLC ( J , 1 ) 

X{ J,2)=XULD( J,2) 

X ( J , 3 ) = XO L C ( J , 3 ) 

207 J=J+1 

GO TO 92 
102 GO TO 1G1 

t NO OF PKCGRAN 
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SUBROUTINE PCTTERCX) 

THIS SUBROUTINE SOLVES A FIVE ELEMENT BANDWIDTH DIAGONAL MATRIX 
WHERE EACH ELEMENT IS ITSELF A 3X3 MATRIX. THE SOLUTION IS BASED 
ON POTTERS METHOD. 

CCMMON/BL1/ ROA, N, PHIO, SOA, DS, DR 

COMMON/ BL5/ A ( 202,3,3) ,B(202,3,3> ,C(202,3,3),D(202,3,3),F(202,?,3), 
1 S ( 202 » 3 ) 

COMMON/BL7/XCLD(2C2,3) , EPS2, EPS 
COMMON/ BL8/ LAM, F10, E20, NU12, NU2I, F T A 
COMMON/ BL9/ IDYN, KM AX , NONLIN, CHAR 
COMMON /BL10/ R ( 1 Cl ) , PHI I 1 Cl ) , DPHI ( 1 01 ) 

DIMENSION XI 202,3) 

DIMENSION AA (3*3) , 68(3,3), CC(3,3I, 00(3,3), EE(3,3), FF(3,3) 

DIMENSION S S ( 3, ) , Fl(3,3), F3(3,3), P(3,3), 0(3,3), AR(3), G(3,3) 

L=2*N 

L1=L-1 

DO 46 M= 1 , 3 

CO 46 MM = 1 ,3 

F 1 ( M , MM ) = 0 • 

46 F3(M,MM)=0. 

IF (NONLIN .EQ. C) GO TO 11 
IF ( C( 1 , 2 , 2 ) .EO. 1.) 0(1, 2,3)= D ( 1 , 2 , 3 ) + X ( 1 , 3 ) 

IF (0(1, 2, 2) .EO. 1.) C( 1,2,1 )= C( 1,2,1 )+X (2, 3 ) 

IF ( B( L , 2, 2 ) .EC. 1.) B(L1,2,1)=B(L1,?,I)+X(L,3) 

IF ( B( L « 2 , 2 ) .FO. 1.) C(L,2, 3 ) = C(L, 2,3) +X(L1 , 1 ) 

DO 23 J=2 , L 1 

C(J,l,3) = C(J,l,3)+(X(J<-2, 3)+X(J,3) )/(4.*ETA) 

E(J,1,3)=E( J,1,?) + (X< J+2, 3)+X( J,3> l/(4.*FTA) 

J = J + 1 
I=( J+l) / 2 
DFT = DPH I ( I ) 

0 SR= COS (PHI ( I ) )/R( I I 
D(J,1,3)=D(J,1 » 3 ) + DFI/(2.*ETA)*X(J,1) 

C(J,l,l)=C(J,l,l)+ DFI /( 2.*ETA) *X( J+l ,3 ) 

R(J,1,3)=B(J,1,3)+ DFI/(?.*ETA)*X(J-2,1 ) 

A(J,1,1)=A(J,1,1)+ OF I / ( 2 .*ETA)*X( J-l ,3) 

C(J,2»3)=D(J,2»? )+ (CSR/( 2.*FTA) +1 ./( ETA*OS) ) *X( J, 1 ) 

C ( J , 2 ,1 ) =C ( J , 2 , 1 ) ♦ (CSR/( 2. + ETA) +1 ./ ( ETA*CS) ) *X( J+! , 3 > 

P(J,2,3) = B(J,2,3)+ (CSR/(2.*ETA)-1./(FTA*DS))*X(J-R,1 ) 

23 A ( J,2 ,1 ) =A( J ,2,1 )+ (CSR/(2.*ETA)-l./(ETA*nS))*X(J-I,3) 

11 CONTINUE 

CO 76 K = 1 , 3 , 2 
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IF (C(l »K»K) «FC • 1.) GO TC 76 
00 7 7 M = I » 3 
P(K,M)= CClfK,n 
C(K,M)= C(IfKfM) 

ClItKfM)- B ( 2 , K » M ) 
r(l,K,M|= C ( 2 t K « M ) 

F(I f K » M ) = 0(2, K,m 
r 1 { K,M) = EC2,K,K) 

P<2, K.M)= P(K 
C ( ? f K *M ) = Q { K t M ) 
r ( ? f K » M ) = o. 

77 b ( 2 ? K »M ) = 0. 

AR { K ) =S ( 1,K) 

S ( 1 , K ) = S ( 2 , K ) 

F(?,K)= AR(K) 

76 CONTINUE 

IF (C(l , 2,2) .FC. 1 • ) GO TC 79 
TO 78 M = 1 * 3 

A < 3 f 3 f M ) 
r (1 f 2 fM)- B( 3 ,3 ,P) 

F(l, ?,M)= C{?,3 f M) 

F I ( 2 , M ) = D f ? , 3 , M > 

A(3,3 f M)= 0. 

P ( 3 ? 3 , M ) = A ( A f 3 » ^ ) 

C(3,3,M)= 

r ( 3 1 3 »m ) = c(4 f 3,n 

F(3 f 3 ,M»= 0(6 ,3 t M 
F 3 ( 3 f M ) = E( A,3,M) 

MA,* f M> = C ( 2 ? 2 t ^ ) 

P(A,3,M)= 0(?,2 f M 
CCA, 3,M)= F ( 2t 2,M 
CCA,3 ,M)= 0. 

F (A,3 ,M)= 0. 

F(2f2»M)— P(2,M) 

C(Z,2 f M)= Q(2»V) 

0(2,2tM)= 0. 

73 FC?,2,M)= 0. 

AR { 2 ) =S ( 1 * 2 ) 

S (!,?)= S ( 3 , 2 ) 

S(3,3>= S ( A , 3 ) 

S(A,3)= S ( 2 * 2 ) 
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S (2,2 > = ARC?) 

79 CONTINUE 

CO 51 M = 1 ,3 
CO 5 2. MM= 1,3 
CCCM, MM)= C ( I , M , MM ) 
00 < M , MM ) =— f) ! I , N f NM) 
52 FF(M f MM l=-F( 1 , M,MM) 
51 SS(M)= S ( 1 , M ) 


CALL 

MA TR I X ( 10, 

3, 

3, 

0, 

CC, 

3, 

CFF» 




CALL 

MATR I X ( 20, 

3, 

3, 

3 f 

CC, 

3, 

DO, 

3, 

00 , 

3) 

CALL 

MATRIX! 20, 

3, 

3, 

3, 

CC, 

3 * 

FF, 

3, 

EE, 

?) 

CALL 

MATRIX! 20, 

3, 

3, 

3, 

CC, 

3, 

FI , 

3, 

FI , 

3) 

CALL 

MATRIX(20, 

3, 

3, 

I t 

CC, 

3, 

SS, 

3, 

SS, 

2 ) 


CO 54 M = 1 , 3 
ro 55 MM=1,3 
D(l f M,MM)= 0D(M,MM) 

F(l f M f MM|= FE<M,MM) 

CC(M,MM)= C(2,M f MN) 

55 R B ( M , MM ) = B(2, M , MM ) 

54 S ( 1 , M )= SS ( M ) 

CALL MA TR I X ( 2 0, 3, 3, 1 , pe, 3, SS, 3, SS, ?) 

CALL MATRIX( 20, 3, 3, 3, BR, 3, DO, 3, G, 3) 

CALL MATRIX! 21, 3, 3, 0, G, 3, CC, 3, G, ?) 
CALL MATRIX(10 f 3, 3, 0, C, 3, GFE) 

CALL MATP I X ( 20, 3, 3, RR, 3, E r , 3, DC, 3) 

CALL MATRIX(20, 3, 3, 3, RB, 3, FI, 3, AA, *) 

CO 58 M = 1 , 3 
00 5^ M M = 1 ,3 

DD( M f MM ) = -CD(M f MM)-D ( 2 , M. , MM! ) 

59 F E ( M , MM ) = AA ( M , M M ) - £ ( 2,M, w M) 

53 SS( M) = -SS(N)+S!2,M> 


CALL 

MATRIX! 20, 

3, % 

3 , 

G, 

3, 

DO , 

3, 

0 0, 

3) 

CALL 

MATRIX! 20, 

3, 3, 

3, 

G, 

3, 

FF, 

3, 

EF, 


CALL 

MATRIX! 20, 

3, 

I , 

G, 

3 , 

SS, 

3, 

SS, 

3) 


LI=L-1 

CO 64 J =3 , L 1 
on 65 M= 1,3 
DO 66 MM=1 , 3 
0( J-l , M, MV)^nn( M f M4) 
E(J-1, M , MV) = FF(M t Mv.) 
A A ( w , MM ) = A ( J , M , M M ) 

B B ( M f M M ) = B!J,M,MM) 
P(M,MM)= D( J-?,M,MM) 
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66 Q!M,MM)= E(J-2,M,MM) 

AR(MI= S ( J-2 , V ) 

65 S! J- 1 * M } = SS(M!) 

CALL MATRIXJ20, 3, 3, 3, AA, 3, P, 3, P, 3) 

CALL MATRIX(20, 3, 3, 3, AA, 3, Q, 3, C, 3) 

CALL MATR IX ( 2 0, 3, 3, I, AA, 3, AR , 3, AR , 3) 


IF 

( J 

• f C • 3) 

CALL 

MATRIX <20* 

3, 

3, 

3, AA, ?, 

FI, 

3, 

AA, 

3) 

IF 

( J 

,FG. 5) 

CALL 

MATRIX! 20, 

3, 

3, 

3, AA, 3, 

F3, 

3 , 

AA, 

3) 

CALL 

MATRIX (21 

, 2, 

3t 0, P, 3 * 

BB 


3, BB, 3* 





IF 

( J 

.EQ. 4) 

CALL 

MATRIX ( 20 , 

3, 

3 * 

3, BB, 3, 

F 3 » 

3, 

AA, 

3) 


CALL MAT R IX ( 2 0, 3, 3, * , BP, 3, DO, 3, DC, 3) 
CALL MATRIX! 20, 3, 3, 3, PP, 3, EE, 3, EE, 3) 
CALL MATR I X ( 20 , 3, 3, 1, BB, 3, SS, 3, SS, 3) 
DO 60 M=1 , 3 
00 67 MM=1,3 

C!M,MM)= DC ( M , MM ) + 0 ( M , MM ) +C I J , M ,MM ) 

DD( M,MM) =-EE ( M , M M ) - D(J,M,MM) 

IF (J .EQ. 3) DC ( M , MM) = DD l M, MM) +AA! M,MM ) 
IF (J .EC. 5) C C ( M , MM) = DD ( M , MM ) + AA ( M , MM ) 
EE(M,MM)= -Et J, M, MM) 

IF (J .EQ. A) F F ( M , MM ) ■= - F ( J , M. , MM ) + A A ( « , MM ) 

67 CONTINUE 

60 5 S ( M ) = S( J,M)-AP(M)-SS(M) 


CAL L 

MATRIX! 1 0, 

3, 

3, 

0, 

G, 

3, 

GEE) 




CALL 

W ATRIX! 20, 

3, 

3, 

3 , 

G, 

3, 

DD, 

3, 

00, 

3) 

CALL 

MATRIX! 20, 

3, 

3, 

3, 

G, 

3, 

EF, 

3, 

FE, 

3) 

CALL 

MATRIX(?0, 

3, 


I , 

G, 

3, 

SS, 


SS, 

31 

IF (J 

• FG . ?) CALL 

MATRIX! 20, 

3, 

3 , 

3, 

G, 

?f 


64 CONTINUE 

CO 68 M=l t 3 

rn 67 M M = 1 , 3 

DILI ,M,MM)= DO(M,MM) 

P(M,MM)= 0( L-2 , M , M M ) 

G { M , MM ) =E ( L -2 , M ,MM ) 

PB(M,MM)= B(L,M,MM> 

A A ( M , M M ) = A(L,M,MM) 

63 CC ( M , MM ) = C(L,M,MM) 

A R ( M ) = S ( L- 2 , M ) 

63 S ( L 1 , M)= SS(M) 

CALL MA TR I X ( 2 0 , 3, 3, i, AA, 3, P, 3, F, 3) 

CALL MATR I X ( 2 0, 3, 3, 3, AA, 3, 0, C, 3) 

CALL MATRIX! 20, 3, 3, l, AA, 3, AR, 7, AR , 3) 

CALL MATRIX!?!, 3, 3, 0, P, 3, BB, 3, BP, 3) 
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CALL MATRIX(20, 3, 3, 1, BB, 3, SS, 3, SS, 3) 

CALL MATR I X I 20? 3* 3, 3, BB, 3, 00, 3, DDt 3) 

CO TO M=l,3 
CO 71 M M = 1 , 3 

71 G ( M , M M ) = OD<M t MM*Q(M,MM) +C(L,M f MM) 

70 S S ( M ) = S(L,M)-AFm-SS(M) 

CALL MATRIX! 10, 3, 3, 0, G, 3, GFE) 

CALL MATR I X ( 2 0 , 3, 3, 1, G f 3, SS, 3, SS, 3) 

CO 30 M— 1 f 3 
CO 32 MM=1,3 

32 CD ( M , MM ) = 0 ( L 1 , M , MM ) 

30 X ( L , M ) = SS < M ) 

CALL MATRIX! 20, 3, 3, It DD, 3, SS, 3, SS, 3) 

CO 33 M= 1 , 3 

33 X ( LI »M)= SS!M)+S!L1,M) 

CO 34 K=2tLl 

J = L~K 

DO 35 M = 1 ,3 
CO 3 6 M M = 1 1 3 
F!M,MM)= D ( J , M , N N ) 

36 G(M, MM) = E( J,M,NM) 
tfi ( M ) = X< J + l y M) 

35 SS(M)= X< J + 2 ,M) 

CALL MATRIX(20, 3, 3, 1, P, ?t AR, 3, AR, ?) 

CALL MATRIX(20, 3, 3, 1, G» 3, SS, 3, SS, 3) 

CO 39 M= 1,3 

39 X { J , M ) = A R ( M ) + S S ( M ) + S(JtM) 

IF (J .GT. 3) GC TO 34 

IF (J • FQ • 2) GC TO 34 
CO 38 Nl = l ,3 

33 5 $ ( M ) = X(J43,M) 

IF (J .EQ. 1) CALL MATRIX! 20, 3, 3, 1, FI, 3, SS, 3, SS, 3) 

IF {J .EQ. 3) CALL MATRIX(20, 3, 3, 1, F?, SS, * , SS, 3) 

ro 40 m=i ,3 

40 X ( J , M ) - X ( J , M ) - $ S ( M ) 

34 CONTINUE 

00 22 J= 1 , L 
CO 22 M= 1,3 

22 XOLO! J, M ) = X( J , M ) 4 XOLD ( J , M ) 

RETURN 

E NO OF S UBRGU T I N F PCTTFR 
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SL.JiNli-JT 1NE ABCCrStX, I, NCYCLE) 

C THIS S L Is rs L L T IKE DEFINES ILL ELEMENTS OF THt A, B, C, D AND E 

C NAT R ICES AS WELL AS THE S VECTOR WHICH IS EQUIVALENT TO THE Q 

C V tC TOR I \ THt LSERS MANUAL 

bcMHul/iiL 1/ RCA » N, J h I G* SCA, DS , DR 
C Cj T- v OK /OL 3 / PP , I ACC * IMT 
CJ^'L H/bLA/PPS 

CDiHJU/OLS/ A( ?C2,3,3) , 0 (202, 3,3 ),C ( 202,3, 3 ) , D< 2C2 , 3 , 3 ) , E < 2 02 , 3 , 3 ) , 
1 5(202,3) 

C-JMMuN/ PL 6/ Cll, C 1 2, C22, Kll, K12, K22, Dll, D12, D22 
COM i / Oi\ / r> L 7 /XC LC(?0?,3 ) , EPS2, EPS 
CJM.VUU/r.LB/ LAP, tlO, E20, NU12, NU2 1 , ETA 
CLH,uii/liL9/ ICW, K MAX , NCNL IN, CHAR 
CCTHi . /;< L 10/ P ( 1C1 ),PHI( 101 ) ,OPHI ( 101 ) 

OJMHClN/h LI 1/ ALTAI , ALFA2, Tl, T2, ITEM' 

CiiM-Vi/RL 12/KCLM 
DI.-ILKS U:x X(2C2, 3) 


kEAL 

NLI2, N UP 1 r 

Kilt 

K 12 t 

K22 

KEAL 

NL2, 3 , M2» M3, 

K612 

» K 3 


KL AL 

L Av 




J= 2 

*i I-i ) 




CALL 

STIMlt LAM, fcl), 

f 20, 

NU12 

, NU2 1 ) 


CUK = (Cl 1*0 1 1-K 1 1**2) /LAM**2 
ETANUS= ETA/( 1 .-NU2 1*NU 12 ) 

ETAL= ETA/LAV 

GL2 = CCK*LAM**2 

GL A= GL2*LAR**2 

iM12= (C12*Cll-Kll*K12)/GL2 

M 1 2 = ( K12*C11-K11*C12 )/COK 

t 12= ( 2.*C 12*Kll*Kl2-C12*C12*Dll-Cll*t<12*K12)/GL2 

KB12= ( C 1 2 * K 1 1*C12-C12*C11*K12<-K12*K12*K11-K12*C11*D12)/GL2 

N3= (K1 2*C 1 1-K 1 1*D12 )/GL4 

M3 = ( D 1 2*C 1 L — K 1 1*K 1 2 ) /GL2 

t 3=(K 1 1*K 1 2*K 12- C 12*D 1 1*K 12 *C 12*01 2* K U-C 1 1 *D 12*K 1 2 ) /GL4 

K3= ( 2 . *0 1 <:*K 11 *K 12- ) 1 1*K 1 2*K 12-C 1 L*U 12*0 12 ) /GL A 

UE= ( 0 1 1*C 12- K 1 1* K 1 2 ) /GL2 

UK- ( K l 1 * L 12-C11*K12)/GL2 

Bc= (K ll*C 12-C11*K12) /GL2 

BK= (C11*012-K11*K12)/GL2 

R I = P ( I ) 

H I = PHI ( I ) 

0 F I = UP H I ( I ) 
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CS= CCS < FI ) 

5N = SI N ( F T ) 

C SP= C S/ R I 
SKR=SN/RI 

TN1=F10*FTANLS*< ALFA1 + NU1 2*AL F A2 ) * T07(I ) 
3N2 = F20 * ET A KU S* ( A IF A 2 + MU? 1 * AL F AI ) * T07 ( M 
TM1=F10*ETAL *< ALFA1 + \IU12*ALFA2 l*TZ07(I > 
TM2=F20*FTAl. * ( AL F A 2 +MU2 1 *AL F A 1 )*TZD7U ) 
A ( J , 1 , 1 ) =0* 

A ( J , 1 , 2 ) = C • 

A ( J 1 1 ,3 ) =0 • 

MJ,2il)=C. 

A ( J , 2 , 2 ) =0 • 

M J f 2,3>=0. 

A { J , 3 » 1 ) = 0 • 

A ( J , 3 ,2 )=0. 

A ( J , 3 * 3 ) =0. 

E( J, 1,1 ) = — D1 1 / ( GL 2* 2 • ) 

P< J,1,2)=C. 

B{ J, 1 ,3 )=Kll/(CCK*2. ) 

F ( J , 2 ,1 )=0. 

P( Jt’t?)=0. 

F ( J, )= 0. 

P ( J , 3 , 1 )=-Kll/(GL2*?. I 
P ( J , 3 , 2 ) = 0 . 

P( J,3,3 )=C11 /(CCK*2, ) 

C(J,1,1>= -1./DS+IE*CS R/2. 

C ( J , 1 , 2 ) = ( DF I +L E* SNR )/2. 

C(J,!,3h CSR*UK/2. 
r (Jf ?fl )= -DFl/2. 

C(J,2,?)= -l./DS 
C<J,2»31= -.5 
C(J,?,1 )= B F 5 P /? • 

C( J, 3,2)= BF*S!VP/2. 

C ( J , 3 , 3 ) = -l./DS + BK+CSH/?. 
n(J,l,l >=-Dll/<GL?*2. ) 

D ( J , 1,2)= 0. 

C { J, 1 ,3 )= K 1 1 /<CCK*2 . ) 

CU, 2,1 ) =0 • 

0!J, ?,2I = 0. 
n ( Jf 2 ,3 )= 0. 

C(J,3 f l)=-Kll/(GL2*2.) 

C ( J, 3,2 )= 0. 
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D(J,3,3)=C11/(CDK*2.) 

F(J,1,1)= l./DS+UF*CSR/2. 

E(J,1»2)= (DFI+LE*SNR) /2. 

E(J,1,31= CSR *UK / 2 . 

F { J, 2 ,1 )= -DFI/2. 

F<J,2,2)= l./DS 
E ( J » 2 » 3 ) = -.5 
F(Jt3*l>- BE*C$R/2. 

F(J,3,2)= BE*SNR/2 • 

E < J * 3 *3 ) = l./BS+BK*CSR/2. 

S(J»l)=011/(l AM**2*CDK)*TN1 -K11/C0K*TM1 
$ ( J » 2 ) = 0 • 

S(Jt 3)= Kll/ <LAM**?*CDK)*TN1-C11 /C0K*TM1 
IF (NONLIN .EQ. CJ GO TO 10 

S(J,1)= IX< J + 2,3) +XU,3))**2/ETA/8. 

10 J=J + 1 

A ( J , 1 , 1 ) = -l./DS+CSR*(l.-N12)/2. 

A ( J * 1,2)= DFI/2. 

A(J,1,3) =— CSR* M 12/2. 

A ( J * 2 ,1 ) = -IDFI+SNP* N 1 2 ) / 2 . 

A(J*2»2)= (-1./CS+CSR/2.) 

A ( J , 2 , 3 )=-SNR*M12/2. 

A(J,3,1)=— CSP*N2/2. 

A(J,3,2)= ~.5/LAM>**2 
A (J ,3 ,3) =-l. /DS + CSR*( 1 .-M3I/2. 

P(J,1,1)= -( E12+CR2 )*CSR*CSR/2. 

P(J,1,*)= -<E12+C2?l*CSR*SNR/2. 

P ( J » 1 * 3 1 = (KB12+K?2)*CSR*CSP/2. 

P. (J, 2,1 )= -( E12+C22 ) *0 SR* SNR/2. 

F(J,?,2>= -( E12-»C22 )*SNR*SNR/2. 

P ( J , 2 ,3 ) = (K312+K2?)*C SR* SNR/2. 

PI J, 3,1 )= -CSR*CSR*(E3*K?2/LAN!**2> / 2. 

0 I J, ? ,2 > = -C SR* SNR* < F 3*K? 2 /LAM** 2 ) / 2 . 

P ( J, 3,3 )= CSR*CSR*IK3+D??/LAM**2 ) /2. 
C(J,1.,1> = l./DS*CSP*(l.-N12)/2. 

C(J,1,?)= DFI/2. 

C I J , 1 ,3 ) =— CSR* M12/2. 

C(J,?,1)= -(OFI+SNR* N12)/2. 

C ( J , 2 , 2 I = ( 1./CS+CSR/2.) 

Cl J,2,3)=-SNR* M 12/’. 

C ( J , 3 , 1 ) =-CSP *N 3 / 2 . 

C ( J , 3 ,? ) = -.5/LAM**? 

C(J»3»3>= I./OS+CSP*Il.-M3)/2. 
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U(J,1,1)= -( E12 + C22)*CSR*CSR/2 . 

D(J* 1*21= — (E12+C22)*CSR*SNR/2. 

D(J, 1,3) = ( KB 12+K22) *CSR*CSR/2 . 

U < J « 2 1 1 )= -( E12+C22) *CSR*SNR/2. 

U( J »2 » 2 )= -( E12+C22)*SNR*SNR/2. 

0( J«2, 3)- (K612+K22)*CSR*SNR/2. 

UI J * 3, l ) = -CSR*CSR*(E3+K22/LAM**2) /2. 

0( J,3,2)= — CS R*SNR* (E3 + K22/LAM<=*2 )/2. 

D(J, 3,3)= CSR*CSR*(K3+D22/LAM**2)/2. 

E(J,1,1)=0. 

E(J,1,2)=C. 

E ( J , 1 , 3 ) =0 . 

E ( J , 2 , 1) = C. 

E (J *2 , 2 ) = C . 

E(J,2, 3)=C« 

E ( J , 3 , 1 ) = C. 

F { J,3,2)=C. 
t( J,3, 3 )=C. 

S(J,1)= -PS(I) +CSR*( N12*TN1+M12*TM1-TN2 > 

S ( J *2 ) = — P L ( I ) +SM^(N12*TN1+M12*TM1— TN2) 

S ( J ,3) = CSF* (N3*TN1 +M3*TM1-TR2 ) 

IP ( I DYR . N 6 . 1 ) GO TO 83 
If ( I A C C . fcC. 1) GO TO 82 
A L= ^./EPS2*I ( I ) 

IF ( KUliNT .EC. INIT) A L=0 . 

IF ( KUUNT .EG. I.NIT + 1) AL= 6 . /£ PS2*T ( I ) 

B( J ,1 , 1 ) = 8 <J, 1,1 ) —A 1/ 2 . 

D< J,l, 1 ) = D ( J, 1,1) — AL/ 2 . 
b(J,2,2)=8(J,2,2)-AL/2. 

D( J,2,2)=D( J,2,2)-AL/2. 

83 CONTINUE 

IP ( 1 Li Y P .LG. 1) CALL CYIMAKIC ( S ,N, J, I, NCYCLF) 
t>2 IP (NJMIN .FC. 0) GO TO 84 

3(J,l)= S(J,1 )-0f 1/FTA*(X(J*1, 3 ) *X (J » 1 ) + X ( J- 1 , 3 ) *X (J-2 , 1) )/2. 
S(J,2)=S(J,2)-CSR/FTA*(X(J+l,3)*X(J,l)+X(J-l,3)*X(J-2,l))/2.- 
L I X ( J + 1 , 3 ) *X ( J , L )-X(J-l,3)*XU-2,l))/ETA/0S 
34 CONTINUE 
RE r Li K 0 

fc NO UP $08 ROUT INt AHC.DPS 



noon 
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SUBROUTINE GFCM7Y (NTYPF) 

THIS SUBROUTINE CALCULATES THE GEOMETRY AT THE STATION MIDPOINT ( 1-1 /2> 
IF NTYPE = A THE USER MUST DEFINE P(I), PH I ( I ) AND DPMI (I) 36 
CORRESPOND WITH HIS DESIRED GECMFTRY. THESE DEFINITIONS SHOULD B p 
l OCA TED BFTWEEN THF LABEL STATEMENTS 50 ANC 60. 

CCMM0N/BL1/ rga, n, phio, soa, ds, dr 
C C M M r N / BL 1 0/ R 1 1 Cl ) ,PUI< 101 ) ,DPHI ( 101 I 
COMMON /BL 1 A/ DS2 , DEL, RPHI, AN 
DATA PI / 3.1415S26535B9793 / 

RPHT = PHIO* PI / 1 P.0.0 
Nl= N-l 
AN =NI 

IF (NTYPF . NE . 1 ) GO TO 10 

C THF SHELL IS A CYLINDFR WHEN NTYPE IS 1 

CS = SOA/ AN 
CO 20 I = 1,N 
P ( I ) = ROA 

PHI (T ) = PPH I 

20 C PH I C I ) = 0.0 

00 TO 60 

10 IF (NTYPE • N E • 2 ) GO TO 30 

C THE SHELL IS A CCNF WIEN NTYPE IS ? 

CS = SOA/ AN 
r R = DS * CCS (RFH I ) 

P ( I ) = RO A-DP / 2 . 

FHI ( 1 ) = PPH I 
CPHim = O.C 
CO 40 I = 2. N 
R(I) = R ( I- 1 ) + CR 

FHI ( I) = RPHI 


40 

CPH I (I ) = 

0.0 



CO TO 60 



30 

IF (NTYPF 

.NE. 

3) GC TO 50 

C 

THE SHELL 

IS A 

SPHERICAL CAP WHEN NTYPE IS 3 


SOA = R OA*RPH I 
CS = SOA/ AN 
OF I =R PH I /AN 
A0R=1 ./ROA 
CP^I ( 1 ) =ACP 
PHI ( 1 ) = -DFI/2. 

F (1) = RO A*S I N ( P H I ( 1 ) ) 

CO 70 1=2, N 

FHI ( I) = PHI ( I -1 ) + DEI 
F ( I ) = ROA * SIN (PHI ( I ) ) 

70 CPU I ( I) =ACR 
GO TO 60 

50 IF (NTYPE .NE. 4) GC TO 60 
60 RETURN 

END OF SUBROUTINE GFCMTY 
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SUBROUTINE ROCCMLRCL, L BCR t X, ITER) 

PQCON DEFINES TEE BOUNDARY CONDITIONS TO PE USED. 

COMMON/ RL 1/ RCA, N, PHIO, SOA, DS, DR 

CCMMON/BL2 /CMEGALI? ,3) , AL AMD AL ( 3 » 3 ) * OMFG AR ( 3 , ? ) , 

1 SL ( 3 ) , SR C 3 > 

COMMON/ BL5/A( 202 ,3,3) , B(202,3 ,3) ,C(202,3,?)»D( 2C2,3i 
1 S I 202 , 3 ) 

COMMON /BL10/ R ( 1 01 ) , P I I ( 1 Cl ) , OP H I ( 1 0 1 ) 

DTMFNSI ON X( ?C2,3) 

DATA PI / 3.1415 9 26535R9793 / 

L = 2*N 
l 1= L-l 

cn l j= l , 3 

CO 2 K=l ,3 
A(l, J,K>=0. 

BIT * J * K I = 0 • 

C ( 1 , J »K ) = C. 

Cll, J,K)=0. 

F ( 1 ,J ,K)=C. 

A(?*N,J ,K)=0. 

P(2*N,J ,K)=0. 

C (2*N,J , K 1 = 0. 
r (2*N,J,K)=0. 

2 F(?*N ,J ,K)=C. 

S(l)Jl-SL(Jl 

1 S(?*N,J)=SR(J) 

IF (L8CL .NF. 1) CO TO 3 
IF ( I TER .EO. 1 ) PRINT 4 

4 FORMA T{ 40X , * TF E LEFT BOUNCARY CONDITION IS A ROLF 
0(1,1, T)=l. 

C ( 1 , 2 , ? ) = 1 . 

0(1,3,31=1. 

CO TO 9 

3 IF (LBCL .NF. 2) CO TO 5 
IF ( ITER .EG. 1 I PRINT 6 

e, FC P M A T( 40X , * TF t LEFT BOUNDARY CONDITION IS PINNt-D* 
D ( 1 , 1 , 1 > = 1 . 

C( 1 , ? ,2 > = 1 . 

C( 1,3,31=1. 

GO TO 9 

5 IF (LBCL .NF. 3) CO TO 7 
IF (ITER .FQ. 1) PRINT 3 


AlAMOAR (3, 
3 ) , E ( 202 , 3 


POINT* ) 
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8 FORMAT* 40X, *TFE LEFT BOUNDARY CCNCITICN IS FIXED*) 

cc 1 , 1 a )=i. 

DUt2»2l*l. 

CI1, 3,31=1. 

GO TO 9 

7 IF { L8CL *NF • 4) GO TO 19 
IF (ITER .FQ. 1) PRINT 10 

10 FORMAT C 40X » *TF£ LEFT BOUNDARY CONDITION IS F R F F * ) 
Ct!il»U=i. 

CC 1,2, 2>“1. 

C C 1 «3 f 3)=l. 

S(lf21= S(l,2!-X(lf U*X(2*3) 

DC TO 9 

19 IF CLBCL •NE. 5) GO TO 9 
IF (ITER .FQ. U PRINT 20 

20 FORMAT* 40X » * TF F LEFT BOUNDARY IS AN ELASTIC CONSTRAINT*! 

DO 21 M=1 * 3 
00 21 MM=1,3 

CCl,M,MM!=GMEGAL*M f MM) 

21 D(i f M,MWJ*Al AMCAl(MrMM) 

9 IF (LBCR .NE. I) GO TO 11 
IF ( ITER .FQ. 1) PRINT 12 

12 FCRMATC40X, *TFE RIGHT BOUNDARY CONDITION IS A PfLF POINT*) 
C*2*N»ltl)— 1* 

F C2*N»2 »2)=1. 

CC2*Nt3 ♦?)=!. 

GO TO 17 

11 IF (LBCR .NE. 2) GO TO 13 
IF (ITER .FQ. 1) PRINT 14 

14 FORMAT* 40X, * THE RIGHT BOUNDARY CONDITION IS PINNED*) 

C ( ?*N 1 1 1 1 ) -1 • 

C ( 2*N » 2 » 2 ) ” I . 
p ( ? * N 1 3 t 3 ) = I . 

GO TO 17 

13 IF (LBCR .NE. 3) CO TO 15 
IF * I TER .EG. 1 5 PRINT 14 

16 FORMAT* 40X, *ThF RIGHT BOUNDARY CONDITION IS FIXED* ) 

C(?*Ntl 1 1)«! . 

C(2*N»2f2)=l. 

C(?*N,3 ,3)=I . 

GO TO 17 

15 IF (LBCR .NE. 4) CO TO 22 
IF (ITER .EG. 1) PRINT 13 
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18 FQRMATC40X, *TFE RIGHT BOUNDARY CONDITION IS FREE-) 

E(2*N,1 ,11=1. 

8f 2*N,2,2)=1. 

E { 2*N ? 3 » 3 )=1 . 

S{L,2)=S(L,2)-X{L1,1)*X{L,3) 

GO TO 17 

22 IF { L BC R .NE. 5) GO TO 17 
IF (ITER .EQ. 1) PRINT 23 

23 FORMAT! 40X, *TFE RIGHT BOUNDARY IS AN ELASTIC CONSTRAINT*) 

DO 24 M=1 1 3 

CO 24 MM=1 » 3 

B{ L, M» MM)=ONEGAF MM) 

24 C(L, M, MM)=AL AMCAR(M,MM) 

17 RETURN 

FND OF BOCON 


APPENDIX F 


SUBROUTINE FUNCT ( X , M 

C THIS SUBROUTINE CALCULATES THE NEGATIVE OF S = AX + P X+GX +0 X+E X-S 

( CMM0N/BL5/A( 20? ,3 , 3 ) , B C 2 02,3 , 3 ) ,C ( 2 0? , 3 , 3 ) , D( 20? , 3 , 3 ) , E ( 202 , 3 , 3 ) , 
I S( 202,3) 

DIMENSION A A ( 3 , 3 ) , BB(3,3>, CC(3,3), 00(3,3), FF(3,3) 

DIMENSION X A ( 3 ) , XR<3), XC(3), XD13), XF(3) 

OINFNSION X ( 202 , 3 ) 

L = ?*N 
l 1=L-1 

no i j=i ,l 

IF ( J .NE. 1 ) or. TO 2 
C 0 3 y= 1 , 3 
XA(M) =0. 

X P ( M ) =0 . 

X C ( M ) = X ( 1,M) 

XD(M )=X ( 2,M) 

3 XF(M)=X(3,M) 

CO TO 10 

2 IF < J . NE. 2 ) GC TO 4 
TO 5 M= 1 * 3 
X A ( M ) =0 , 

>B(M) =X( 1 ,M> 

XC ( M ) = X ( 2 , M ) 

XC <M)=X ( 3 ,M) 

5 XE(M)=X(4,M) 

GO TO 10 

4 IF (J .NE. L) GC TO 6 
no 7 M= 1 , 3 

X A ( M ) =X ( L— 2 , V ) 

X B ( M ) =X( L-1,M) 

XC(M)=X(L,M) 

X 0 ( M ) =0. 

T XE(M)=0. 

GO TO 10 

6 IF (J .NE. ( L - 1 ) ) GC TO 3 
TO 9 M=l,3 

X A ( M ) = X ( L — 3 , M ) 

XE ( H ) =X ( L— 2 , M ) 

XC ( M ) =X ( L— 1 , P ) 

XC(M) =X( L ,M) 

9 X E ( M ) =0 • 

GO TO 10 
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3 CO 11 M=l,3 
>A(M>= X l J-2 » M) 

X B ( M } = X ( J — 1 * M ) 
XC(M)=X( J,M) 
XC(M)=X< 

11 XR{ M)=X ( J+2, M) 

10 00 12 ’1 = 1,3 

CO 12 M M= 1 , 3 
AA(M,MM)=AI 
PB(M,'1M1 = B( ) 

CC(M,MM)=C(J»M,PM) 
C0(M,MMJ=0( 

12 FE«M,MH = E( J,M,NK ) 


CALL 

MATRIX! 20* 

3, 

3, 

1 , 

AA, 

3, 

XA , 

3 , 

XA, 

3 ) 

CALL 

MATRIX! 20, 

3 f 

3, 

I , 

PR, 

3, 

XB , 

3, 

XB , 

31 

CALL 

MATRIX! 20, 

3 , 

3 , 

1 , 

CC, 

3 , 

xc , 

3, 

xc. 

?) 

CALL 

MATR IX< 20, 

3, 

3, 

1 , 

DO, 

3, 

XD , 

3 , 

xc, 

3) 

CALL 

MATRIX(20, 

3, 


1, 

EE, 

3, 

XE» 

3, 

XF, 



CO 1 * M= 1 , 3 

13 S(J,M)= S(J,M) -Xfl(M)- XRIM)-XC(M)- XDIM)-Xfc(M) 

1 CONTINUE 
RETURN 

END OF SUBROUTINE FUNCT 
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SUBROUTINE ANSWER SI NTYPE , X, LBCL, LBCR) 

C THIS SUBROUTINE CALCULATES AND PRINTS N-THETA AND M-THFTA. 

C THE LOAD VECTORS ARE ALSO PRINTED WITH THESE VALUES. 

CCMMON/BL1/ ROA, N * PHIO, SOA , DS, DR 

CCMM0N/BL6/ Cll, C12, C22» Kll, K12, K22, Dll, D12, D22 
COMMON/ BL8/ LAM, E10, E20, NU12, NU21, ETA 
COMMON / BLIO/ R ( 1C1 1 , PH I ( 101 1 , DP HI ( 101 ) 

COMMON/BL11/ALFA1, ALFA2, Tl, T2, ITFMP 
COMMON /BL14/ DS2 , DEL, RPHI, AN 
DIMENSION X( 202 ,3 I 
PEAL LAM, Kll, K12, K22 
PEAL N12, N3 , M12, M3, KB12, K3 
PEAL NU1 2 , NU2 1 
CALL GE C MTV (NTYPE) 

L = 2*N 

1 FORMAT! /5X*ST ATICN*7X*NTHETA*14X*MTHETA*14X *R AD I A L L 0AD*7X*MF R I D I 
IAN LOAD* ) 

PRINT 1 
CO 2 1=1, N 
CALL STIF(I) 

FTAL= ETA/LAM 

f TANU S= ETA/(1.-NL21*NU12 ) 

lNl = E10*ETANliS*(ALFAl+NU12*ALFA2)* TDZ ( I ) 

TN2 = E20*ETANUS*( ALF A2+NU2 1*AL F A1 )* T0Z( I ) 

TM1=E 10*ETAL * ( A LF A 1+NU1 2*AL F A2 ) *T ZOZ ( I ) 

TM2=F?0*ETAL * ( ALF A2+NU2 1 *ALF A1 ) *TZDZ ( I ) 

L I = X ( 2* I ,1) 

H I = X (2*1 ,2) 

PI=X( 2*1 ,3) 

IF( I.NE.ll GO TC 1C 
F I = PH I ( 1 ) +DFL 
DF I = 0PHK1) 

RT= IRC 1 J+R< 211/2. 

CO TO 11 

10 F t = PH I ( I ) +OFL 
C F I = DP H I ( I ) 

PI= R (I ) +DS2 *CC S ( PHI ( I ) ) 

11 CONTINUE 

IF { I . NE. 1 I GC TO 3 
IF (LBCL .NF. 1) GO TO 12 

LP= (-3 ,*X(2 , 1 )+A .*X (4, 1 1 -XI 6, 1) ) /< 2.*DS) 

PP= (-3.*X(2,3)+A.*X(4,3)-X(6,3) )/(2.*DS) 
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GO TO 15 

3 IF (I • NE« N) GC TO 12 
IF ( LBCR .NF. i) GO TO 12 

LP= ( X(L-4 t n-A.*ML-2,l)+3.*X(L,ll )/<2.*DS) 

8P= ( X( L-4 f 3 )-4#*ML-2» 31+3* *XIL » 3 ) )/< 2.*DS) 

15 AK22=-8P 

EM22=UP+DFI*k I+BI**2/(?.*ETA) 

GO TO 13 

12 S NR= SINCFD/RI 
CSR=COS (FI)/PI 

EM22= CSR*UI+SNF *WI 
AK?2 = -CSR*BI 

13 COK = < C11*D11-Kll**2)/LAM**2 
GL 2= CD K*LAM* *2 

C L4= GL 2 *L AM * *2 

M2= <C12*D11-K11*K12> /GL2 

N 12= <K12*C11-K1I*C12 )/CDK 

El 2= <2.*C12*K11*K12-C1 2* C 12*0 1 1 -C 1 l *K 1 2*K 1 2 ) / GL 2 
K B1 2= ( C12*K1 i*D12-C 12*01 1*K 1 2+K 1 2*K 1 2*K 3 1-K1 2 *C 1 1 * D1 2 ) /GL? 

N 3= 1K1?*D11-K 11*012 ) /GL4 
P3= (D12*C11-K11*K1?) /GL2 

E3=(KX1*K12*K12- C12*011*K12 + C 1 2*D 1 >*K 1 1-C 1 1 *D 1 ?*K 2 ? ) /G L4 
K 3= (2. *D12*K11*K12-011*K12*K12-C1 1*D1 2*D12)/GL4 

AN22=N12*X( 2*1-1 » 1 )+M12*X < 2*1-1, 3) +(E12+C?2)*E M 22-M K 3 1 2 + K 2 Z ) * A K 2 ? 
2 + N12*TN1 + M12*TR1-TN2 

A M22= N3 *X(2*I-lt 1) + M3 *X C 2* I - 1 , 3 ) ♦ < F 3 + K2? / L AM** 2 ) * FM 2 2 

1 +(K3+D22/LAM**2)*AK?2 

2 + N3*TNl + N ! 3*T^l-TM2 
PPL = P L ( I , N ) 

PPS=P$( l ,N1 

PRINT 6, I, A N2 2 * AM?2, PPL, P PS 
6 FCRM AT( 4X , 15, <i(5X, F15.8)) 

2 CONTINUE 
RETURN 

FND OF ANSWERS SLPRCUTINE 
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SUBROUTINE ST IF I I ) 

C THIS SUBROUTINE DEFINES Oil. C12, C22t Kilt K!2, K?2, Dll, Q12 

C AND 022 • 

CCMM0N/BL6/ Cll, CI2, C2? , Kll, K12, K22, Cll, DI2, D22 

COMMON/ BL8/ LAM, F10. E20, NU12, NL?1, ETA 

CCMMCN/BL12/KCUNT 

C0MM0N/BL1 5/ HO 

c-’EAL MULC , MULD 

PEAL LAM, NU1 2t NU2I, Kll, K12» K22 
C El 0= El/EO AND E2C= E2/E0 

MULC= TID/I1 .-NL12*NU?1) 

NULD = L AM**2*T< I) **2* MI JLC / 12. 

K 11=0 . 

K 12=0 • 

K 2 2=0 • 

f ll=F 10+MUIC 
012= NU12*F1C*MLLC 
f 2 ? = F 2 0 * M U L C 
[]!= C10*MULD 
012= NU1 2*E10*MLLC 
022= E20*MULD 
RETURN 

F NO OF SUBROUTINE STIF 
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SUIJkUCT 1 N E CYMFIC(S, N, J, I, NCYCLE ) 

THIS SLBKCLT INE SETS ALPHA, BETA, GAMMA AMO DELTA FfJK THfc 
BACKWARD TIME DERIVATIVES AND SAVES THt BACKWARD TIME STATIUN 
SULUTICNS IF L AND W. 

C U M MbN/BL 3/ PP » I ACC, I N IT 
CCMMCN/BL7/XOLC( 202,3) , EPS?, EPS 
COMMON / f'L 1 ? /K CENT 

C UMMUN /PL 1 3? X 1 ( 1 0 1 , 2) , X2( IC1 ,2 ) , X3II01,?) 

DIMENSION S ( ?C? , 3 ) 

IF (KGLNT .NE. I NIT) GO TO 1 
CE=0 . 

GA=0. 

UE= l . 

IF (IP(KCLNT) .EC. 1) GO TO A 
If (NCYCL F .CT . I ) GO TO A 
GO TO A 

1 IE ( KUUM .NF. { IN I T+ 1 ) ) GO TO 2 
IF (NCYCLE .CT. 1 ) GO TC 10 

DU 6 M = 1 » 2 

o XI ( I,y ) = < XCl D ( J + l ,F )+ X LLD ( J— 1 1 M ) )/2. 
lu BE= -6./EP$?*T( 1 ) 

GA = -6. /EPS *T( I ) 

DE = -2 . 

GO TO A 

2 IF (KCLNT . F F . (IMIOII GO TC 3 
if (NCYCLE .CT. 1 ) GJ TO 11 

DO 7 M = t , 2 

X?( I , M )=X 1 ( I ,F ) 

7 XI ( I ,m ) = ( XCLD ( J + l ,F )+ XGED ( J— 1 , M) ) /2 . 

11 b F.= — A./EPS?*T(I ) 

GA= 2 ./fc PS2*T ( I ) 

Ofc-=-l . 

Gu ru A 

5 IF (NCYCLE .GT. 1) CO TC 12 
DO B M= 1,2 
X 3 ( I , M ) ■= X 2 ( 1 , F ) 

X2( I » M ) = X 1 ( I,F) 

o XI ( I,M) = (XCLD(J+1,M)+ XOLD ( J-1.,M) ) /2. 

12 b L= -t>. /LP S2*T ( I ) 

GA= A . / 1 P S 2 * T ( I ) 

Dh= -1 ./EPS?* K I ) 

DO 9 M= 1,2 

9 S ( J , M ) •= S ( J,MtliE*Xl(l , ,v )<-GA*X2( I , M ) + X 3 ( I , M ) *DL 
l< t T Lk INI 

ENO OF DYN AFT C 
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SUBROUTINE INITI4UX, LBCL I 
THIS SUBROUTINE INITIALIZES THE X AND X DOT VECTORS 

THIS SUBROUTINE IS REQUIRFD IF X OR XOOT (INITIAL CONDITIONS! APE 
CTHER THAN ZERC FOR DYNAMIC PROBLEMS. 

COMMON/ BL1/ ROA, N, PHIO, SOA , OS, OR 

C0MM0N/BL6/ Cll, Cl 2 , C22 , Kilt K12, K22t Dll, Dl?, 022 
COMMON/ BL 8/ LAM, E10, E20, NU12, NU21, ETA 
CCMMCN/BL9/ IDYM, KMflX, NCNLIN, CPAR 
COMMON /BL10/ R ( 1C1 I ,PH I { 1C1 1 ,OPHI ( 101 I 
COMMON/ BL13/X 1( 1C1, 2 ) , X2(101,2), X3(101,2) 

DIMENSION X( 202 , 3 ) 

REAL LAM, Kll, K12, K22 
CO 16 1=1 ,N 
J = 2*I 

X(J, 1 )= DV( 1 » I > 

X{J,2)=DV(2,I ) 

X 2 ( 1 , 1) =DV( 3, I) 

16 X 2 1 1 ,2)=DV(4, I) 

no l i=l , n 

J=2*I 

IF (I .NE. 1) GC TO 3 

WP=(-3.*X(2*2)*4.*X(4,2)-X(6,2)!/(?.*0S) 

CC TO 5 

3 IE (I .NE. N) GC TO A 
WP=(X(J-4,2)-4.*X(J-2, ? M-?.*X(J,?»)/(2.*DS> 

CFI =DPH I ( M + ( OPE T (M! -OPHT (N-1H/2. 

GO TO 1 

4 RP=(X<J+2,2J-X( J-2,2) )/(2.*0S) 

5 rEI= ( DPHI ( n + DPFI ( 1 + 11 ) /? . 

1 X ( J , 3 )= WP-DF I*X < J , 1 I 

DO 2 T = 1 , K 
J = ?*I 

CALL STIFU, LAM, F10, E10, NU12, NU21! 

IF (I .NE. 1) GC. TO 6 

LP=(-3. *X(2 , 1 ) *4 .*X ( 4, 1 )-X<6, 1 ) ) /( 2 .*0S) 

PP=( - 7 . *X (2,3‘» + 4 .*X ( 4, 3 )-X (6,31 ) /( ?.*0S) 

GO TO 3 

6 IF (I . NE. M CC TO 7 
UP=(X(J-4,ll-4.*X(J-2,l)+?.*X( J,l> )/(2.*CS) 
PP=(X(J-4,3!-4.*X(J-2,3 M-3.*X( J,3I l/(2.*CS) 

PI=R (N) +CCS (PHI (M )*CS/2« 

F I = PHI(N)+(PHT (N1-PHI ( \|- L > ) / 2. 
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CFI=DPHI (N) + { DPHI (N) -OPHI (N-in/2. 

GO TO 9 

7 LP=< X(J + 2,l )-X( J-2,1 ) ) /<2 .*DS ) 

BP={X(J+2,3)-X( J-2,3) J/(2.*DS) 

3 F I = * R ( I )+R(I + I) )/2. 

F 1= ( PHI ( I l+PHK 1 + 1) ) /2 . 

GFI=( OPHim+OPF KI + IH/?. 

9 IF (LBCL .NE. 1 ) GO TO 10 
IF (I . EQ. 1) GG TO 11 

10 SNP=SIN(F I)/F I 
CSR*C0S ( FI) /R I 
FM22=CSR*X< J , 1 ) + SNR*X(J ,2 ) 

AK22=-CSR*X( J ,3) 

11 FM11=UP+DFI*X( J,2) 

IF (NONLIN .EC. 1) EM11= EMI 1 +X ( J , 3 ) ** 2 / ( 2 . *E T A ) 

AK11=— BP 

IF ILRCL .NE. 1) GO TO 12 
IF (I . NE. 1 ) GC TO 12 
FM2 2 = EM 1 1 
AK22=AK1 1 

12 M = J— 1 

X ( M * 1 >=Cll*EMll+C12*E >i '22 + Kll'f ! AK3 1+K12*AK22 
X ( M , 3 ) = ( K11+EM1 1+K1?*EM2?+011*AK1 1+D12*AK22 )/LAV**? 

IF (LBCL .NF. 1) GO TO 17 
IF (I .EQ. 1) X ( M , 3 ) =0 . 

GO TO 2 

17 X(M»2)=CSR*(X(M|?)-(K12*EM11 + K22*E‘'122 + D12*AK’2 + P?P j «'AK22)) 

2 CONTINUE 
GO 13 1 = 1, N 
J=2*I-1 

IF (I .NE. 1) GC TO 14 

X(J,?)=XIJ,2)+LAM**2*<-3.*X(2t3M-4.*X(4,?)-X(6,3))/(’.*QSI 
CO TO 13 

14 IF (I .NE. N) GC TO 15 

X(J,2)=X(J,?)+L AM** 2*13. *X(J t 3>-4.*X(J-2,3)+X(J- 4, ?))/(’. 
ro TO 13 

15 X ( J » 2 )= X( J ,2 ) +LAM**?*( X( J+2, 3 >-X( J-2, ?) ) /( 2.*DS> 

13 GCNTINUE 
RETURN 

END OF INITIAL 
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f UNCTION PL( I ) 

COMMON/ RL 1/ POA, N, PHIO* SQA , D$, DR 
rCMMC-N/BL3/PP f IACC, IN IT 

COMMON/ B L 8/ LAM, FiO, F20, NU12, NU21 t CIA 

CCMMON/BLI2/KGUM 

PEAL LAM 

PEAL NU12 f NU21 

FC.L = 2.*LAN*E7 A/( 3 .* ( 1 . -MJ 1 2 *NU2 I ) ) **•?*< T( I )/POA )**2*F 10 
PLL=PP*PCL 
F L- PL L 

TF { KCUNT • F G . C ) PL = D. 

IF (KO'JNT .GT.2C ) PL = 0. 

IF (TACC .NF. 1 ) GC IN 1 

I<= ( K PUNT • FQ • C) PL= PLL 
IF (KCUNT . FC .2C ) PL= 0. 

1 CONTINUE 
RETURN 
f NO 


FLNCT IPN PS( I) 

CCM MCN/BL1/ R 0 A , N, PHIO f SOA , OS, OR 
CC’ v Mf N/SL3/PP , IACC, TMIT 
CCMM0N/BL4/PP S 

C OMMC N/ B L P / LAM, CIO, =20, NU12, NU?L, ETA 

CCNMCN7BL12/KCLNT 

F S=PP S 

f FTUPN 

END 


FUNCTION IP (KCUNT) 

IP=0 

TF ( K TU NT * F 0 • 0) IP=! 

IF (KOUNT • F 0 • 20 IP = 1 
RETURN 
LND 
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FUNCTION DV(Ntl) 


D V ( 1 * 

I ) 

SFTS 

T F F 

INITIAL 

U 

DI SPLACEPFNT 

D V ( 2 » 

I ) 

SETS 

TEE 

INITIAL 

w 

DISPLACFMFNT 

ovn, 

T ) 

SETS 

TEE 

INITIAL 

u 

VFLCCI TY 

DV (A, 

I ) 

SETS 

TEE 

INITIAL 

vi 

VELQC. I TY 


c v=o* 

RETURN 
END CF DV 


FUNCTION T ( I ) 

CCMM0N/BL15/ HC 

C TCI) MUST BE N CNP I N c N S I CNA L THICKNESS, THFREF OR F DIVIDE PY - i 0 
T =1 • / HO 
PFTUPN 
END 


FUNCTION TOZIII 
CCMMON/BLl/ P 0 A » N, 

puro, sca, 

ns, 

DR 

COPMCN/BLl 1/ ALF AT , 

ALFA?, Tl, 

T ?. , 

I TFMP 

I T = 1 TEMP+1 





A N — I T EM P + 1 
H= T ( I) /?. 

TCZ-T1 *TU)+T2*(H**IT-(-H>**IT)/A\ 
RETURN 
F ND 


FUNCT ION TZDZ ( I ) 

COMMON/ 0 L I / RCA, N , PHIO, SOA , OS, DR 
COMMON / BL 1 1 / AL F A 1 , ALFA 9 , Tl, T2, ITEMP 
I T=I TEMP+2 
AN=ITFMP+2 
H=T ( I) /?. 

I ZD Z= T, 2 * ( F** I T-(-H)**m / AN 

RFTURN 

END 
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SAMPLE PRINTOUT 

The printout for the sample problem is as follows: 


SGI VEN 



NTYPE 

= 

3* 

N 

= 

26, 

RO 


0.1E+C3, 

SO 

= 

0.0, 

HO 

= 

0.1E+C1, 

EO 


0.1E+C8, 

El 

= 

0.1E+C6, 

E2 

= 

o.iE+ce, 

NU1 2 


0.3E-»00, 

NU2 1 


0.3E+00, 

SIGO 

- 

0 • IE +01 • 

NONLIN 

= 

1 , 

CONV 


0.1E-C2, 

PHIO 

= 

0.13 RCS4E+02 , 

LBCL 

= 

l, 

LBCR 

= 

3, 

PP 

= 

-0 *6E ♦ OO, 

PPS 

= 

o 

o 

SL 

» 

0.0, 0.0, 0.0 

SR 

- 

o 

o 

o 

o 

o 

o 

CHAR 

- 

0.1E+C3, 

1DYM 

= 

It 

KMAX 

= 

40, 

DTAU 


0.25E+00, 



ri = o.o, 

T2 = 0.0, 

ALFA1 = 0.0, 

ALFA2 = 0.0, 

ITEMP = 0, 

IFREQ = 4, 

ISTART = 0, 

SEND 

SAMPLE PROBLEM FOR A CLAMPEC SFEERICAL CAP WITH LAMBDA - = 


*0 
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—3 

CD 


IF IDYM=0 THE SHELL IS LOADED STATICALLY 
IF ICYM = 1 THE SHELL IS LOADED DYNAMICALLY 
IF I D YM = 2 THE SHELi STATIC BUCKLING LOAD IS CALCULATED 


El/EO a 
NU 1 2 = 
EO/SG * 


FCR THIS RUN, ICYM= 

NUMBER OF STATIUNS = 

MFP IDIAN/REFFRENCE RATIO = 
RADIUS/REFERENCE RATIO = 
THICK/ REF RAD RATIO = 
1.000000 
.300000 
1CCC0000. 000000 


l 

26 

.275926 
1.000000 
.010000 
E2/E0 = 

NU2 1 = 

REF DIST * 


IF NONLIN = 0 UNLY LINEAR TERMS ARE USED 
IF NGNLIN = 1 NLNLINfAR TERMS USED 


FOR THIS RUN NCNL IN = l 


1.000000 

.300000 

100.000000 
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THE GEOMETRY Of TEE SHELL FOLLOWS 


STATION NO 

MERIDIAN CIS 

RADIAL DIS 

ANGLE (RAD) 

CUPVATUPEtDFI/DS) 

1 

0 . 

-5. 6C2033E-C 8 

- 2 .775558 E- 17 

l.OOOOOOE+OO 

2 

1 • 10 37C4E-02 

l. 1C3693E-02 

1.103704E-0? 

1 • OOCOOOE+OO 

3 

2.207409E-C2 

2.207257F-02 

2.2C7409F-02 

! • OOOOOOF + CO 

4 

3. nil 13E-02 

3* 31055 3E-02 

3.31111 3E-02 

l.OOOOOOF+OO 

5 

4.414817E-02 

4.413445E-02 

4. 4 14817E- C2 

l.OOOOOOF+OO 

6 

5* 51 8522E-02 

5. 515799E-02 

5.5 L 8522 E-02 

1. OOOOOOE+OO 

. 7 

6 .622226E-02 

6. 6174 8 2E- 02 

6.622226F-C? 

l.OOOOOOF+OO 

8 

7.725930E-C2 

7.718359E-02 

7 • 7 25°30E-02 

i : oooooqe+oo 

9 

8.829635E-02 

8. 818295E-C2 

8.829635 E-02 

l.OOOOOOF+OO 

10 

9.933339E-02 

9 • 9 171 57E-02 

9 . 9 33 " 39F-02 

1. OOOOOOE+OO 

11 

1 • 1037C4E-01 

1. 101481F-01 

1.1C3704F-01 

l.OOOOOOE+OO 

12 

1.214075E-C1 

1.21U12E-CI 

1. 2 1 4075E- Cl 

l.OOOOOOE+OO 

13 

1 • 724445E- 01 

1 • 320596E-01 

1.324445E-C1 

l.OOOOOOE+OO 

14 

1.43481 6E- 01 

1.429919F-C1 

1.434316F-01 

1 • 000000 F +00 

15 

1 • 54 5 1 86E-01 

1. 539067E-01 

1.5451 8 A F-01 

l.OOOOOOE+OO 

16 

1 • 655 556E- 01 

1.648C29E-01 

1.655556E-G1 

l.OOOOOOE+OO 

IT 

1 .76 5927E-0 1 

I. 756789E-01 

1.765 927F- Cl 

l.OOOOOOF+OO 

18 

1.876297E-01 

1.865335E-01 

1.P76297F-01 

l.OOOOOOE+OO 

19 

1.98 6668E- 0 1 

1 . 97 36 5 5F-C1 

1.986668E-01 

l.OOOOOOE+OO 

20 

2.C97038E-01 

2 • 0 81 73 3E-0 1 

2.097038F-01 

1 .000000 c +00 

21 

2.207409E-C1 

2. 1 99559E-0L 

2. 207^0° E-01 

1 • OOCOOOE+OO 

22 

2*31 7779F- 01 

2.297117E-01 

2.31 777°F- Cl 

l.OOOOOOE+OO 

23 

2.428150E-01 

2. 404396E-01 

2.428150E-01 

1 • OOOOCOE+OO 

24 

2 . 53 8520E- C 1 

2.5 1 13 8 IE— Cl 

2 . 5 38520E-0Z 

I .OOOOOOE+OO 

25 

2.648890E-01 

2.61 806 1E-0 1 

2.64B890F-01 

l.OOOOOOF+OO 

26 

2 . 75926 IF— C 1 

2. 7 24422F-01 

2.759261F-01 

1 .OOOOOOE+OO 


THE MAXIMUM SHELL RISE IS .03782669 

TEE SHELL THICKNESS = 1.00000000 

THE LEFT 5CUNCARY CONDITION IS A POLE PCINT 
THE RIGHT BOUNDARY CONDITION IS FIXED 



THE 

LEFT OMEGA MATRIX 




THE 

LEFT LAMBDA MATRIX 





0 . 


0 . 

0 . 


+ 

i .oooooE+no 

0 . 

0 . 

(UJ 


0 . 

0 . 


1 • OOOCOE+OO 

0 . 

(Cl 

+ 

0 . 

0 . 

0 . 

<wi 

= 

0 . 

0 . 


0 . 

0 . 

(M) 

+ 

0 . 

0 . 

1 . OPOOOE +00 

<P) 

= 

0 . 


THE 

RIGHT OMEGA MATRIX 




THE 

RIGHT L AM 80 A MATRIX 





0 . 


0 . 

0 . 

«M 

♦ 

1 • ODOOOE+OO 

0 . 

0 . 

(U) 

= 

0 . 

0 . 


0 . 

0 . 

(C) 

+ 

0 . 

3 .OOCOOE+OO 

0 . 

( W I 

= 

0 . 

0 . 


0 . 

0 . 

m 

+ 

0 . 

0 . 

1 • CCOOOE +00 

IP) 

= 

0 . 


-3 

CO 


APPENDIX G 



00 

o 


STATION NO 


N-S RESULTANT 


SHEAR FORCE 


M-S RESULTANT 


1 

0. 


2 

0. 


3 

0. 


A 

0. 


5 

0. 


6 

0. 


7 

0. 


8 

0. 


9 

0. 


10 

0. 


11 

0. 


12 

0. 


13 

0. 


14 

0. 


15 

0. 


16 

0. 


17 

0. 


18 

0. 


19 

0. 


20 

0* 


21 

0. 


22 

0. 


23 

0. 


24 

0. 


25 

0. 


26 

0. 


STATION 

NTHETA 


1 

0. 

0. 

2 

0, 

0. 

3 

0. 

0. 

4 

0. 

0. 

5 

0. 

0. 

6 

0. 

0. 

7 

0. 

0. 

8 

0. 

0. 

9 

0. 

0. 

10 

0, 

0. 

11 

0. 

0. 

12 

0. 

0. 

13 

0. 

0. 

14 

0. 

0. 

15 

0. 

0. 

16 

0. 

0. 

17 

0. 

0. 

18 

0. 

0. 

19 

0. 

0. 

20 

0. 

0. 

21 

0. 

c. 

22 

0. 

0. 

23 

0. 

0. 

24 

0. 

0. 

25 

0. 

0* 

26 

0. 

0. 


0. C. 

0 . 0 . 

0. C. 

0* 0. 

0. C. 

0. C. 

0. 0. 

0, c. 

0. c. 

0. c. 

0 . 0 . 

0 . 0 . 

0 . 0 . 

0 . 0 . 

0 . c. 

0 . 0 . 

0 • c • 

0 . 0 . 

0 . 0 . 

0 . 0 . 

0 . 0 . 

0 . c. 

0 . 0 . 

0 . 0 . 

0. c. 

0 . 0 . 

MTHETA RADIAL LOAD 


0 . 

0 . 

0 . 

0 . 

0 . 

0 . 

0. 

0 . 

0. 

0. 

0. 

0. 

0 . 

0 . 

0 . 

0 . 

0 . 

0 . 

0. 

0 . 

0 . 

0 . 

0 . 

0 . 

0. 

0 . 


ooocooocoooooooooooo 


U-DEFORMAT I ON 


W-DEFOP NATION 


BETA ROTATION 


0. 

0. 

0. 

0. 

0. 

0. 

MERIDIAN 

0. 

0. 

0. 

0. 

0. 

0, 

0. 

0. 

0. 

0. 

0 . 

0 . 

0 . 

0 . 

0. 

0 . 

0 . 

0 . 

0 . 

0 * 

0 . 

0. 

0 * 

0 . 

0 . 

0* 


0. 

0. 

0. 

0. 

0. 

0. 

0. 

0. 

0, 

0. 

0. 

0. 

0. 

0, 

0. 

0. 

0. 

0. 

0. 

0. 

0. 

0. 

0. 

0. 

0. 

0. 

LOAD 


0 . 

0. 

0 . 

0 . 

0 . 

0 . 

0. 

0 . 

0 . 

0. 

0. 

0. 

0 . 

0 . 

0 . 

0. 

0 . 

0 . 

0, 

0. 

0 . 

0. 

0 . 

c. 

0 . 

0 . 
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STATION 

NO 


U DOT 


W COT 

U DOTDOT 

W DOTDOT 

1 


C. 


0. 


0. 

0. 

2 


c. 


C. 


0. 

-7.?fc27303°F+04 

3 


c. 


C. 


C. 

- 7 .?6273039E+04 

4 


c. 


C. 


0. 

-7 , 26273039E+04 

5 


c. 


C. 


c. 

-7. 26273039E+04 

6 


c. 


0. 


c. 

-7. 262 7303 9P + 04 

7 


c. 


0. 


c. 

-7, 26 273039F+04 

S 


c. 


C. 


c. 

-7. ?62 73 0? C T + 04 

9 


c. 


0. 


0. 

-7.26?73039E+04 

10 


c. 


C. 


c. 

-7. 26 27 3039 F +04 

11 


c. 


0. 


0. 

-7. 262 7 3 03 °F + C4 

12 


c. 


0. 


0. 

-7. 262 7303 c E+0^ 

13 


c. 


0. 


c. 

-7.2627303QF+04 

14 


c. 


0. 


0. 

-7. ?627303 QfI + 04 

15 


c. 


c. 


c. 

-7, 26273039E+04 

16 


c. 


0. 


0. 

-7. 26 2 73 03 9F + 04 

17 


c. 


c. 


0. 

- 7 .262730?9 i: + 0 4 

18 


c. 


0. 


c. 

-7. 26273039F+04 

19 


c. 


0. 


0. 

-7. 26 27303 C F +04 

20 


c. 


c. 


c. 

-7.26273039F+04 

21 


c. 
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STATION' NO 


N-S RESULTANT 


SHFAR FORCE 


M-S RESULTANT 


1 

-2. 95360 8CE + C4 

2 

-2.9536507E+C4 

3 

-2.9547457E + C4 

4 

-2.9569002E4C4 

5 

— 2.9604770F+C4 

6 

-2.9659037F+C4 

7 

-2.9735975E4C4 

8 

-2. 983868 6E + C4 

9 

— 2.9968121E+C4 

10 

-3. 012200 1E + C4 

11 

-3. 029391 3E + C4 

12 

-3.0472829F4C4 

13 

-3.0643271E+C4 

14 

-3.0786302E+C4 

15 

-3.0881279E+C4 

16 

-3.09C3078E+C4 

17 

-3.0849245F+C4 

18 

- 3. 0691 51 7E + C4 

19 

-3.04265C4F+C4 

20 

-3.00508C8E+C4 

21 

-2.956636 IE + C4 

22 

-2.8981643E+C4 

23 

-2. 83 137 8 5E + C4 

24 

-2.7 590 54 OE + C 4 

25 

-2. 6350791F + C4 

26 

-2. 614311 5F4C4 


0. 2 • 2678 369F + 04 

1. 0P23427E+C2 2 . 86 c 102 OF + 04 

2* 0C92919E+C? 4 .28 57752F+0 4 

2. 62847 C2F+C? 6 . 3 P 57 895F +04 

2.7904913E+02 P .794^5^ RF +04 

2. 362 5C43E + C? 1 . 100?9<=?F + Q5 

1. 24942 26F+02 1 . 239 ? 3 5RF +C 5 

-5. 7701113F+01 1 .229 10780 +C 5 

-3.0459378E+C2 1 . C05P51 5f +05 

-5, 9597947 E + 02 5 . 16967 59F +04 

-S.9818565F+C2 . 57 5322 8F+04 

-1. 16 49886E + C3 - 1 * 299 7R32E +0 5 
-1. 3424 182 E+03 -2.5491 796F +05 

-1. 3771512E+C3 -3.855271 8F+0 5 

-1. 22757 C8E+03 -5 .C75026 78 +05 

-8.75321 11E+C2 -6.C017726F+05 

-3. 34166 54E + C2 -fc . 439501 ?E +05 
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STATION 

NTHETA 

MTFETA 

RADIAL LOAD 


MERIDIAN LOAD 
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-7, 2627303 9F+ 04 

0. 


2 
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APPENDIX H 

CONVERSION OF U.S. CUSTOMARY UNITS TO SI UNITS 

The International System of Units (SI) was adopted by the Eleventh General Con- 
ference on Weights and Measures in 1960. (See ref. 18.) Conversion factors for the 
units used in this report are given in the following table: 


Physical quantity 

U.S. Customary Unit 

Conversion factor 
(*) 

SI Unit 
(**) 

Length 

Modulus of axial 

in. 

2.54 x 10“ 2 

meter (m) 

stress, elasticity 

psi 

6.895 x 10 3 

newton/meter 2 (N/m2) 

Temperature 

degree Fahrenheit 

K = (°F + 459. 67) /l. 8 

kelvin (K) 


Multiply value given in U.S. Customary Unit by conversion factor to obtain 


equivalent value in SI unit. 

**The prefix giga (G) is used to indicate 10^ units. 
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TABLE 1.- SUBROUTINE DESCRIPTIONS 


FORTRAN name 

Description 

GE0MTY 

Calculates shell geometry, that is, r, <p , and <p * 

B0C0N 

Sets proper boundary conditions 

ABODES 

Calculates coefficients of A, B, C, D, E, and q matrices 

FUNCT 

Calculates negative of fk(zi,Zi_i,s^ in equation (27) 

P0TTER 

Calculates P, Q, and R matrices and then solves for the Newton - 
Raphson corrections 6z 

DYNAMIC 

Calculates a, /3, y, and 6 and stores X m ^_j, X m ^_ 2, 
and X m> £_3 

ANSWERS 

Calculates 1122 and m22 and prints n22, m 22> P> and Ps 

STIF 

Calculates constitutive coefficients in equations (14) 

INITIAL 

Initializes y^ q and y^ q vectors when the initial conditions are 
other than zero 

MATRIX 

Perform basic matrix operations as detailed in reference 16 
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TABLE 2.- GLOSSARY OF FORTRAN VARIABLE NAMES 


Variable 

Program name 

Description 

a 

CHAR 

Reference length 

C ll> C 12’ C 22 

Cll, C12, C22 

Constants in equations (10) and (11) 

Du, D 12> d 22 

Dll, D12, D22 

Constants in equations (12) and (13) 

E 0 

E0 

Reference modulus of elasticity 

Ei/Eq, E 2 /Ej 

E10, E20 

Nondimensional moduli of elasticity 

e ll> e 22 

EMU, EM22 

Membrane strains, used in ANSWERS 

H 0 

H0 

Reference shell thickness 

h 

T(I) 

Nondimensional shell thickness 

3 

ITEMP 

Temperature exponent in equation (17) 

K ll> k 12’ k 22 

Kll, K12, K22 

Constants in equations (10) to (13) 

m 

K0UNT 

Step in time 

n 

N 

Number of stations 

P> Ps 

PL, PS 

Nondimensional loads, see FUNCTION PL and 
FUNCTION PS 

r i 

R(I) 

Nondimensional radial distance at station i— 1/2 
defined in GE0MTY 

s 

SMER 

Nondimensional meridional distance 

x (1st element) 

X(J,1) 

nji when j odd 

x (2d element) 

X(J,2) 

q when j odd 

x (3d element) 

X(J,3) 

mu w hen j odd 

y (1st element) 

X(J,1) 

u when j even 

y (2d element) 

X(J,2) 

w when j even 

y (3d element) 

X(J ,3) 

(3 when j even 

z 

XOLD 

Matrix of unknown 

Sz 

X 

Matrix of Newton-Raphson corrections to z 

«1» «2 

ALFA1, ALFA2 

Constants in equations (16) 

a, % y, 6 

AL, BE, GA, DE 

Constants in equation (21) 
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TABLE 2.- GLOSSARY OF FORTRAN VARIABLE NAMES - Concluded 


Variable 

Program name 

Description 

V 

ETA 

E 0 /a 

K ll> K 22 

AK11, AK22 

Nondimensional principal curvatures 

"12’ "21 

NU12, NU21 

Poisson's ratios 

X 

LAM 

Ratio of H 0/ / a 

*-s 

LAMS 

Shell parameter 

,m .m 

h ’ l 2 

TM1, TM2 

Thermal moment resultants 

t n t n 

l l> l 2 

TNI, TN2 

Thermal force resultants 

At 

ATAU 

Nondimensional time increment 

As 

DS 

Nondimensional meridional increment 

<Pi 

PHI (I) 

Colatitude angle at i- 1/2 


DPHI(I) 

Colatitude angle change at i- 1/2 
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(a) Membrane and trans- 
verse force resultants. 


0 



(b) Moment resultants. 



Figure 2.- Positive sense of forces, moments, loads, 
displacements, and rotations. 
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(a) Shell meridian divided into (b) Location of i- 1/2 station 

equal increments. 

Figure 3.- Locations of boundaries, stations, and midpoints of shell and 

shell meridian increments. 



Figure 4.- Spherical cap geometry. 



= 5 

Results in reference 12 / 



Dimensionless time, t 


Figure 5.- Deflection response of clamped spherical cap for various 
uniformly distributed pressures. (Load duration is from t = 0 
to r = 5.) 



Maximum inward deflection, 



0 .2 .4 .6 .8 

Load parameter, P* 


Figure 6.- Dynamic buckling of clamped spherical cap subjected to 
uniformly distributed step pulses. (Load duration is from 
t = 0 to r = 5.) 


Axial stress at outside surface, psi 



Axial stress, Gl\l/m 





Circumferential stress at outside surface, psi 



Circumferential stress, GN/m 2 
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